1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Utility functions for setting up SHOCKTUBE
6
7 #include "../qfunctions/shocktube.h"
8
9 #include <ceed.h>
10 #include <petscdm.h>
11
12 #include <navierstokes.h>
13
PRINT_SHOCKTUBE(Honee honee,ProblemData problem,AppCtx app_ctx)14 static PetscErrorCode PRINT_SHOCKTUBE(Honee honee, ProblemData problem, AppCtx app_ctx) {
15 MPI_Comm comm = honee->comm;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscPrintf(comm,
19 " Problem:\n"
20 " Problem Name : %s\n",
21 app_ctx->problem_name));
22 PetscFunctionReturn(PETSC_SUCCESS);
23 }
24
25 static const char *const component_names[] = {"Density", "MomentumX", "MomentumY", "MomentumZ", "TotalEnergy"};
26
NS_SHOCKTUBE(ProblemData problem,DM dm,void * ctx)27 PetscErrorCode NS_SHOCKTUBE(ProblemData problem, DM dm, void *ctx) {
28 SetupContextShock setup_context;
29 Honee honee = *(Honee *)ctx;
30 MPI_Comm comm = honee->comm;
31 Ceed ceed = honee->ceed;
32 PetscBool implicit;
33 PetscBool yzb;
34 PetscInt stab;
35 ShockTubeContext shocktube_ctx;
36 CeedQFunctionContext shocktube_qfctx, ics_qfctx;
37
38 PetscFunctionBeginUser;
39 PetscCall(PetscNew(&setup_context));
40 PetscCall(PetscNew(&shocktube_ctx));
41
42 // ------------------------------------------------------
43 // Create the QFunction context
44 // ------------------------------------------------------
45 // Driver section initial conditions
46 CeedScalar P_high = 1.0; // Pa
47 CeedScalar rho_high = 1.0; // kg/m^3
48 // Driven section initial conditions
49 CeedScalar P_low = 0.1; // Pa
50 CeedScalar rho_low = 0.125; // kg/m^3
51 // Stabilization parameter
52 CeedScalar c_tau = 0.5; // -, based on Hughes et al (2010)
53 // Tuning parameters for the YZB shock capturing
54 CeedScalar Cyzb = 0.1; // -, used in approximation of (Na),x
55 CeedScalar Byzb = 2.0; // -, 1 for smooth shocks
56 // 2 for sharp shocks
57 PetscReal domain_min[3], domain_max[3], domain_size[3];
58 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
59 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
60
61 // ------------------------------------------------------
62 // Command line Options
63 // ------------------------------------------------------
64 PetscOptionsBegin(comm, NULL, "Options for SHOCKTUBE problem", NULL);
65
66 // -- Numerical formulation options
67 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL));
68 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL));
69 PetscCall(PetscOptionsScalar("-c_tau", "Stabilization constant", NULL, c_tau, &c_tau, NULL));
70 PetscCall(PetscOptionsBool("-yzb", "Use YZB discontinuity capturing", NULL, yzb = PETSC_FALSE, &yzb, NULL));
71
72 // -- Warnings
73 if (stab == STAB_SUPG) {
74 PetscCall(PetscPrintf(comm, "Warning! -stab supg not implemented for the shocktube problem. \n"));
75 }
76 if (yzb && implicit) {
77 PetscCall(PetscPrintf(comm, "Warning! -yzb only implemented for explicit timestepping. \n"));
78 }
79
80 PetscOptionsEnd();
81
82 // ------------------------------------------------------
83 // Set up the QFunction context
84 // ------------------------------------------------------
85 // -- Scale variables to desired units
86 CeedScalar mid_point = 0.5 * (domain_size[0] + domain_min[0]);
87
88 // -- Setup Context
89 setup_context->mid_point = mid_point;
90 setup_context->time = 0.0;
91 setup_context->P_high = P_high;
92 setup_context->rho_high = rho_high;
93 setup_context->P_low = P_low;
94 setup_context->rho_low = rho_low;
95
96 // -- QFunction Context
97 honee->phys->implicit = implicit;
98 shocktube_ctx->implicit = implicit;
99 shocktube_ctx->stabilization = stab;
100 shocktube_ctx->yzb = yzb;
101 shocktube_ctx->Cyzb = Cyzb;
102 shocktube_ctx->Byzb = Byzb;
103 shocktube_ctx->c_tau = c_tau;
104
105 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &ics_qfctx));
106 PetscCallCeed(ceed, CeedQFunctionContextSetData(ics_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*setup_context), setup_context));
107 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(ics_qfctx, CEED_MEM_HOST, FreeContextPetsc));
108
109 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &shocktube_qfctx));
110 PetscCallCeed(ceed, CeedQFunctionContextSetData(shocktube_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*shocktube_ctx), shocktube_ctx));
111 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(shocktube_qfctx, CEED_MEM_HOST, FreeContextPetsc));
112
113 // ------------------------------------------------------
114 // SET UP SHOCKTUBE
115 // ------------------------------------------------------
116 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsShockTube, .qf_loc = ICsShockTube_loc, .qfctx = ics_qfctx};
117 problem->apply_vol_rhs = (HoneeQFSpec){.qf_func_ptr = EulerShockTube, .qf_loc = EulerShockTube_loc, .qfctx = shocktube_qfctx};
118 problem->num_comps_jac_data = 0;
119 problem->compute_exact_solution_error = PETSC_FALSE;
120 problem->print_info = PRINT_SHOCKTUBE;
121
122 problem->num_components = 5;
123 PetscCall(PetscMalloc1(problem->num_components, &problem->component_names));
124 for (PetscInt i = 0; i < 5; i++) PetscCall(PetscStrallocpy(component_names[i], &problem->component_names[i]));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127