1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Utility functions for setting up SHOCKTUBE
10
11 #include "../qfunctions/shocktube.h"
12
13 #include <ceed.h>
14 #include <petscdm.h>
15
16 #include "../navierstokes.h"
17
NS_SHOCKTUBE(ProblemData problem,DM dm,void * ctx,SimpleBC bc)18 PetscErrorCode NS_SHOCKTUBE(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
19 SetupContextShock setup_context;
20 User user = *(User *)ctx;
21 MPI_Comm comm = user->comm;
22 Ceed ceed = user->ceed;
23 PetscBool implicit;
24 PetscBool yzb;
25 PetscInt stab;
26 ShockTubeContext shocktube_ctx;
27 CeedQFunctionContext shocktube_context;
28
29 PetscFunctionBeginUser;
30 PetscCall(PetscCalloc1(1, &setup_context));
31 PetscCall(PetscCalloc1(1, &shocktube_ctx));
32
33 // ------------------------------------------------------
34 // SET UP SHOCKTUBE
35 // ------------------------------------------------------
36 problem->dim = 3;
37 problem->ics.qfunction = ICsShockTube;
38 problem->ics.qfunction_loc = ICsShockTube_loc;
39 problem->apply_vol_rhs.qfunction = EulerShockTube;
40 problem->apply_vol_rhs.qfunction_loc = EulerShockTube_loc;
41 problem->apply_vol_ifunction.qfunction = NULL;
42 problem->apply_vol_ifunction.qfunction_loc = NULL;
43 problem->compute_exact_solution_error = PETSC_FALSE;
44 problem->print_info = PRINT_SHOCKTUBE;
45
46 // ------------------------------------------------------
47 // Create the libCEED context
48 // ------------------------------------------------------
49 // Driver section initial conditions
50 CeedScalar P_high = 1.0; // Pa
51 CeedScalar rho_high = 1.0; // kg/m^3
52 // Driven section initial conditions
53 CeedScalar P_low = 0.1; // Pa
54 CeedScalar rho_low = 0.125; // kg/m^3
55 // Stabilization parameter
56 CeedScalar c_tau = 0.5; // -, based on Hughes et al (2010)
57 // Tuning parameters for the YZB shock capturing
58 CeedScalar Cyzb = 0.1; // -, used in approximation of (Na),x
59 CeedScalar Byzb = 2.0; // -, 1 for smooth shocks
60 // 2 for sharp shocks
61 PetscReal domain_min[3], domain_max[3], domain_size[3];
62 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
63 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
64
65 // ------------------------------------------------------
66 // Create the PETSc context
67 // ------------------------------------------------------
68 PetscScalar meter = 1e-2; // 1 meter in scaled length units
69 PetscScalar second = 1e-2; // 1 second in scaled time units
70
71 // ------------------------------------------------------
72 // Command line Options
73 // ------------------------------------------------------
74 PetscOptionsBegin(comm, NULL, "Options for SHOCKTUBE problem", NULL);
75
76 // -- Numerical formulation options
77 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL));
78 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL));
79 PetscCall(PetscOptionsScalar("-c_tau", "Stabilization constant", NULL, c_tau, &c_tau, NULL));
80 PetscCall(PetscOptionsBool("-yzb", "Use YZB discontinuity capturing", NULL, yzb = PETSC_FALSE, &yzb, NULL));
81
82 // -- Units
83 PetscCall(PetscOptionsScalar("-units_meter", "1 meter in scaled length units", NULL, meter, &meter, NULL));
84 meter = fabs(meter);
85 PetscCall(PetscOptionsScalar("-units_second", "1 second in scaled time units", NULL, second, &second, NULL));
86 second = fabs(second);
87
88 // -- Warnings
89 if (stab == STAB_SUPG) {
90 PetscCall(PetscPrintf(comm, "Warning! -stab supg not implemented for the shocktube problem. \n"));
91 }
92 if (yzb && implicit) {
93 PetscCall(PetscPrintf(comm, "Warning! -yzb only implemented for explicit timestepping. \n"));
94 }
95
96 PetscOptionsEnd();
97
98 // ------------------------------------------------------
99 // Set up the PETSc context
100 // ------------------------------------------------------
101 user->units->meter = meter;
102 user->units->second = second;
103
104 // ------------------------------------------------------
105 // Set up the libCEED context
106 // ------------------------------------------------------
107 // -- Scale variables to desired units
108 for (PetscInt i = 0; i < 3; i++) {
109 domain_size[i] *= meter;
110 domain_min[i] *= meter;
111 }
112 problem->dm_scale = meter;
113 CeedScalar mid_point = 0.5 * (domain_size[0] + domain_min[0]);
114
115 // -- Setup Context
116 setup_context->mid_point = mid_point;
117 setup_context->time = 0.0;
118 setup_context->P_high = P_high;
119 setup_context->rho_high = rho_high;
120 setup_context->P_low = P_low;
121 setup_context->rho_low = rho_low;
122
123 // -- QFunction Context
124 user->phys->implicit = implicit;
125 shocktube_ctx->implicit = implicit;
126 shocktube_ctx->stabilization = stab;
127 shocktube_ctx->yzb = yzb;
128 shocktube_ctx->Cyzb = Cyzb;
129 shocktube_ctx->Byzb = Byzb;
130 shocktube_ctx->c_tau = c_tau;
131
132 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &problem->ics.qfunction_context));
133 PetscCallCeed(ceed,
134 CeedQFunctionContextSetData(problem->ics.qfunction_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*setup_context), setup_context));
135 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(problem->ics.qfunction_context, CEED_MEM_HOST, FreeContextPetsc));
136
137 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &shocktube_context));
138 PetscCallCeed(ceed, CeedQFunctionContextSetData(shocktube_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*shocktube_ctx), shocktube_ctx));
139 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(shocktube_context, CEED_MEM_HOST, FreeContextPetsc));
140 problem->apply_vol_rhs.qfunction_context = shocktube_context;
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
PRINT_SHOCKTUBE(User user,ProblemData problem,AppCtx app_ctx)144 PetscErrorCode PRINT_SHOCKTUBE(User user, ProblemData problem, AppCtx app_ctx) {
145 MPI_Comm comm = user->comm;
146
147 PetscFunctionBeginUser;
148 PetscCall(PetscPrintf(comm,
149 " Problem:\n"
150 " Problem Name : %s\n",
151 app_ctx->problem_name));
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154