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 EULER_VORTEX
10
11 #include "../qfunctions/eulervortex.h"
12
13 #include <ceed.h>
14 #include <petscdm.h>
15
16 #include "../navierstokes.h"
17
NS_EULER_VORTEX(ProblemData problem,DM dm,void * ctx,SimpleBC bc)18 PetscErrorCode NS_EULER_VORTEX(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
19 EulerTestType euler_test;
20 User user = *(User *)ctx;
21 StabilizationType stab;
22 MPI_Comm comm = user->comm;
23 Ceed ceed = user->ceed;
24 PetscBool implicit;
25 EulerContext euler_ctx;
26 CeedQFunctionContext euler_context;
27
28 PetscFunctionBeginUser;
29 PetscCall(PetscCalloc1(1, &euler_ctx));
30
31 // ------------------------------------------------------
32 // SET UP DENSITY_CURRENT
33 // ------------------------------------------------------
34 problem->dim = 3;
35 problem->ics.qfunction = ICsEuler;
36 problem->ics.qfunction_loc = ICsEuler_loc;
37 problem->apply_vol_rhs.qfunction = Euler;
38 problem->apply_vol_rhs.qfunction_loc = Euler_loc;
39 problem->apply_vol_ifunction.qfunction = IFunction_Euler;
40 problem->apply_vol_ifunction.qfunction_loc = IFunction_Euler_loc;
41 problem->apply_inflow.qfunction = TravelingVortex_Inflow;
42 problem->apply_inflow.qfunction_loc = TravelingVortex_Inflow_loc;
43 problem->apply_outflow.qfunction = Euler_Outflow;
44 problem->apply_outflow.qfunction_loc = Euler_Outflow_loc;
45 problem->compute_exact_solution_error = PETSC_TRUE;
46 problem->print_info = PRINT_EULER_VORTEX;
47
48 // ------------------------------------------------------
49 // Create the libCEED context
50 // ------------------------------------------------------
51 CeedScalar vortex_strength = 5.; // -
52 CeedScalar c_tau = 0.5; // -
53 // c_tau = 0.5 is reported as "optimal" in Hughes et al 2010
54 PetscReal center[3], // m
55 mean_velocity[3] = {1., 1., 0}; // m/s
56 PetscReal domain_min[3], domain_max[3], domain_size[3];
57 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
58 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
59
60 // ------------------------------------------------------
61 // Create the PETSc context
62 // ------------------------------------------------------
63 PetscScalar meter = 1e-2; // 1 meter in scaled length units
64 PetscScalar second = 1e-2; // 1 second in scaled time units
65
66 // ------------------------------------------------------
67 // Command line Options
68 // ------------------------------------------------------
69 PetscOptionsBegin(comm, NULL, "Options for EULER_VORTEX problem", NULL);
70 // -- Physics
71 PetscCall(PetscOptionsScalar("-vortex_strength", "Strength of Vortex", NULL, vortex_strength, &vortex_strength, NULL));
72 PetscInt n = problem->dim;
73 PetscBool user_velocity;
74 PetscCall(PetscOptionsRealArray("-mean_velocity", "Background velocity vector", NULL, mean_velocity, &n, &user_velocity));
75 for (PetscInt i = 0; i < 3; i++) center[i] = .5 * domain_size[i];
76 n = problem->dim;
77 PetscCall(PetscOptionsRealArray("-center", "Location of vortex center", NULL, center, &n, NULL));
78 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL));
79 PetscCall(PetscOptionsEnum("-euler_test", "Euler test option", NULL, EulerTestTypes, (PetscEnum)(euler_test = EULER_TEST_ISENTROPIC_VORTEX),
80 (PetscEnum *)&euler_test, NULL));
81 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL));
82 PetscCall(PetscOptionsScalar("-c_tau", "Stabilization constant", NULL, c_tau, &c_tau, NULL));
83 // -- Units
84 PetscCall(PetscOptionsScalar("-units_meter", "1 meter in scaled length units", NULL, meter, &meter, NULL));
85 meter = fabs(meter);
86 PetscCall(PetscOptionsScalar("-units_second", "1 second in scaled time units", NULL, second, &second, NULL));
87 second = fabs(second);
88
89 // -- Warnings
90 if (stab == STAB_SUPG && !implicit) {
91 PetscCall(PetscPrintf(comm, "Warning! Use -stab supg only with -implicit\n"));
92 }
93 if (user_velocity && (euler_test == EULER_TEST_1 || euler_test == EULER_TEST_3)) {
94 PetscCall(PetscPrintf(comm, "Warning! Background velocity vector for -euler_test t1 and -euler_test t3 is (0,0,0)\n"));
95 }
96
97 PetscOptionsEnd();
98
99 // ------------------------------------------------------
100 // Set up the PETSc context
101 // ------------------------------------------------------
102 user->units->meter = meter;
103 user->units->second = second;
104
105 // ------------------------------------------------------
106 // Set up the libCEED context
107 // ------------------------------------------------------
108 // -- Scale variables to desired units
109 for (PetscInt i = 0; i < 3; i++) {
110 center[i] *= meter;
111 domain_size[i] *= meter;
112 mean_velocity[i] *= (meter / second);
113 }
114 problem->dm_scale = meter;
115
116 // -- QFunction Context
117 user->phys->implicit = implicit;
118 euler_ctx->curr_time = 0.;
119 euler_ctx->implicit = implicit;
120 euler_ctx->euler_test = euler_test;
121 euler_ctx->center[0] = center[0];
122 euler_ctx->center[1] = center[1];
123 euler_ctx->center[2] = center[2];
124 euler_ctx->vortex_strength = vortex_strength;
125 euler_ctx->c_tau = c_tau;
126 euler_ctx->mean_velocity[0] = mean_velocity[0];
127 euler_ctx->mean_velocity[1] = mean_velocity[1];
128 euler_ctx->mean_velocity[2] = mean_velocity[2];
129 euler_ctx->stabilization = stab;
130
131 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &euler_context));
132 PetscCallCeed(ceed, CeedQFunctionContextSetData(euler_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*euler_ctx), euler_ctx));
133 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(euler_context, CEED_MEM_HOST, FreeContextPetsc));
134 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(euler_context, "solution time", offsetof(struct EulerContext_, curr_time), 1,
135 "Physical time of the solution"));
136 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(euler_context, &problem->ics.qfunction_context));
137 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(euler_context, &problem->apply_vol_rhs.qfunction_context));
138 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(euler_context, &problem->apply_vol_ifunction.qfunction_context));
139 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(euler_context, &problem->apply_inflow.qfunction_context));
140 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(euler_context, &problem->apply_outflow.qfunction_context));
141 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&euler_context));
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
PRINT_EULER_VORTEX(User user,ProblemData problem,AppCtx app_ctx)145 PetscErrorCode PRINT_EULER_VORTEX(User user, ProblemData problem, AppCtx app_ctx) {
146 MPI_Comm comm = user->comm;
147 Ceed ceed = user->ceed;
148 EulerContext euler_ctx;
149
150 PetscFunctionBeginUser;
151 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->ics.qfunction_context, CEED_MEM_HOST, &euler_ctx));
152 PetscCall(PetscPrintf(comm,
153 " Problem:\n"
154 " Problem Name : %s\n"
155 " Test Case : %s\n"
156 " Background Velocity : %f,%f,%f\n"
157 " Stabilization : %s\n",
158 app_ctx->problem_name, EulerTestTypes[euler_ctx->euler_test], euler_ctx->mean_velocity[0], euler_ctx->mean_velocity[1],
159 euler_ctx->mean_velocity[2], StabilizationTypes[euler_ctx->stabilization]));
160
161 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->ics.qfunction_context, &euler_ctx));
162 PetscFunctionReturn(PETSC_SUCCESS);
163 }
164