xref: /libCEED/examples/fluids/problems/eulervortex.c (revision e6225c4737d2b5358b9c1c5b7c13d86915d842f6)
1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3 // reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 /// @file
18 /// Utility functions for setting up EULER_VORTEX
19 
20 #include "../navierstokes.h"
21 #include "../qfunctions/setupgeo.h"
22 #include "../qfunctions/eulervortex.h"
23 
24 PetscErrorCode NS_EULER_VORTEX(ProblemData *problem, void *setup_ctx,
25                                void *ctx) {
26   EulerTestType     euler_test;
27   SetupContext      setup_context = *(SetupContext *)setup_ctx;
28   User              user = *(User *)ctx;
29   StabilizationType stab;
30   MPI_Comm          comm = PETSC_COMM_WORLD;
31   PetscBool         implicit;
32   PetscBool         has_curr_time = PETSC_TRUE;
33   PetscBool         has_neumann = PETSC_TRUE;
34   PetscInt          ierr;
35   PetscFunctionBeginUser;
36 
37   ierr = PetscCalloc1(1, &user->phys->euler_ctx); CHKERRQ(ierr);
38 
39   // ------------------------------------------------------
40   //               SET UP DENSITY_CURRENT
41   // ------------------------------------------------------
42   problem->dim                     = 3;
43   problem->q_data_size_vol         = 10;
44   problem->q_data_size_sur         = 4;
45   problem->setup_vol               = Setup;
46   problem->setup_vol_loc           = Setup_loc;
47   problem->setup_sur               = SetupBoundary;
48   problem->setup_sur_loc           = SetupBoundary_loc;
49   problem->ics                     = ICsEuler;
50   problem->ics_loc                 = ICsEuler_loc;
51   problem->apply_vol_rhs           = Euler;
52   problem->apply_vol_rhs_loc       = Euler_loc;
53   problem->apply_vol_ifunction     = IFunction_Euler;
54   problem->apply_vol_ifunction_loc = IFunction_Euler_loc;
55   problem->apply_sur               = Euler_Sur;
56   problem->apply_sur_loc           = Euler_Sur_loc;
57   problem->bc                      = Exact_Euler;
58   problem->setup_ctx               = SetupContext_EULER_VORTEX;
59   problem->bc_func                 = BC_EULER_VORTEX;
60   problem->non_zero_time           = PETSC_TRUE;
61   problem->print_info              = PRINT_EULER_VORTEX;
62 
63   // ------------------------------------------------------
64   //             Create the libCEED context
65   // ------------------------------------------------------
66   CeedScalar vortex_strength = 5.;    // -
67   PetscScalar lx             = 1000.; // m
68   PetscScalar ly             = 1000.; // m
69   PetscScalar lz             = 1.;    // m
70   PetscReal center[3], mean_velocity[3] = {1., 1., 0};
71 
72   // ------------------------------------------------------
73   //             Create the PETSc context
74   // ------------------------------------------------------
75   PetscScalar meter    = 1e-2; // 1 meter in scaled length units
76   PetscScalar second   = 1e-2; // 1 second in scaled time units
77 
78   // ------------------------------------------------------
79   //              Command line Options
80   // ------------------------------------------------------
81   ierr = PetscOptionsBegin(comm, NULL, "Options for EULER_VORTEX problem",
82                            NULL); CHKERRQ(ierr);
83   // -- Physics
84   ierr = PetscOptionsScalar("-vortex_strength", "Strength of Vortex",
85                             NULL, vortex_strength, &vortex_strength, NULL);
86   CHKERRQ(ierr);
87   PetscInt n = problem->dim;
88   PetscBool user_velocity;
89   ierr = PetscOptionsRealArray("-mean_velocity", "Background velocity vector",
90                                NULL, mean_velocity, &n, &user_velocity);
91   CHKERRQ(ierr);
92   ierr = PetscOptionsScalar("-lx", "Length scale in x direction",
93                             NULL, lx, &lx, NULL); CHKERRQ(ierr);
94   ierr = PetscOptionsScalar("-ly", "Length scale in y direction",
95                             NULL, ly, &ly, NULL); CHKERRQ(ierr);
96   ierr = PetscOptionsScalar("-lz", "Length scale in z direction",
97                             NULL, lz, &lz, NULL); CHKERRQ(ierr);
98   n = problem->dim;
99   center[0] = 0.5 * lx;
100   center[1] = 0.5 * ly;
101   center[2] = 0.5 * lz;
102   ierr = PetscOptionsRealArray("-center", "Location of vortex center",
103                                NULL, center, &n, NULL); CHKERRQ(ierr);
104   ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation",
105                           NULL, implicit=PETSC_FALSE, &implicit, NULL);
106   CHKERRQ(ierr);
107   ierr = PetscOptionsEnum("-euler_test", "Euler test option", NULL,
108                           EulerTestTypes, (PetscEnum)(euler_test = EULER_TEST_ISENTROPIC_VORTEX),
109                           (PetscEnum *)&euler_test, NULL); CHKERRQ(ierr);
110   ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL,
111                           StabilizationTypes, (PetscEnum)(stab = STAB_NONE),
112                           (PetscEnum *)&stab, NULL); CHKERRQ(ierr);
113   // -- Units
114   ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units",
115                             NULL, meter, &meter, NULL); CHKERRQ(ierr);
116   meter = fabs(meter);
117   ierr = PetscOptionsScalar("-units_second","1 second in scaled time units",
118                             NULL, second, &second, NULL); CHKERRQ(ierr);
119   second = fabs(second);
120 
121   // -- Warnings
122   if (stab == STAB_SUPG && !implicit) {
123     ierr = PetscPrintf(comm,
124                        "Warning! Use -stab supg only with -implicit\n");
125     CHKERRQ(ierr);
126   }
127   if (user_velocity && (euler_test == EULER_TEST_1
128                         || euler_test == EULER_TEST_3)) {
129     ierr = PetscPrintf(comm,
130                        "Warning! Background velocity vector for -euler_test t1 and -euler_test t3 is (0,0,0)\n");
131     CHKERRQ(ierr);
132   }
133 
134   ierr = PetscOptionsEnd(); CHKERRQ(ierr);
135 
136   // ------------------------------------------------------
137   //           Set up the PETSc context
138   // ------------------------------------------------------
139   user->units->meter  = meter;
140   user->units->second = second;
141 
142   // ------------------------------------------------------
143   //           Set up the libCEED context
144   // ------------------------------------------------------
145   // -- Scale variables to desired units
146   lx = fabs(lx) * meter;
147   ly = fabs(ly) * meter;
148   lz = fabs(lz) * meter;
149   for (int i=0; i<3; i++) center[i] *= meter;
150 
151   // -- Setup Context
152   setup_context->lx        = lx;
153   setup_context->ly        = ly;
154   setup_context->lz        = lz;
155   setup_context->center[0] = center[0];
156   setup_context->center[1] = center[1];
157   setup_context->center[2] = center[2];
158   setup_context->time      = 0;
159 
160   // -- QFunction Context
161   user->phys->stab                        = stab;
162   user->phys->euler_test                  = euler_test;
163   user->phys->implicit                    = implicit;
164   user->phys->has_curr_time               = has_curr_time;
165   user->phys->has_neumann                 = has_neumann;
166   user->phys->euler_ctx->curr_time        = 0.;
167   user->phys->euler_ctx->implicit         = implicit;
168   user->phys->euler_ctx->euler_test       = euler_test;
169   user->phys->euler_ctx->center[0]        = center[0];
170   user->phys->euler_ctx->center[1]        = center[1];
171   user->phys->euler_ctx->center[2]        = center[2];
172   user->phys->euler_ctx->vortex_strength  = vortex_strength;
173   user->phys->euler_ctx->mean_velocity[0] = mean_velocity[0];
174   user->phys->euler_ctx->mean_velocity[1] = mean_velocity[1];
175   user->phys->euler_ctx->mean_velocity[2] = mean_velocity[2];
176   user->phys->euler_ctx->stabilization    = stab;
177 
178   PetscFunctionReturn(0);
179 }
180 
181 PetscErrorCode SetupContext_EULER_VORTEX(Ceed ceed, CeedData ceed_data,
182     AppCtx app_ctx, SetupContext setup_ctx, Physics phys) {
183   PetscFunctionBeginUser;
184 
185   CeedQFunctionContextCreate(ceed, &ceed_data->setup_context);
186   CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST,
187                               CEED_USE_POINTER,
188                               sizeof(*setup_ctx), setup_ctx);
189   CeedQFunctionContextCreate(ceed, &ceed_data->euler_context);
190   CeedQFunctionContextSetData(ceed_data->euler_context, CEED_MEM_HOST,
191                               CEED_USE_POINTER,
192                               sizeof(*phys->euler_ctx), phys->euler_ctx);
193   if (ceed_data->qf_ics)
194     CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->euler_context);
195   if (ceed_data->qf_apply_sur)
196     CeedQFunctionSetContext(ceed_data->qf_apply_sur, ceed_data->euler_context);
197   if (ceed_data->qf_rhs_vol)
198     CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->euler_context);
199   if (ceed_data->qf_ifunction_vol)
200     CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, ceed_data->euler_context);
201 
202   PetscFunctionReturn(0);
203 }
204 
205 PetscErrorCode BC_EULER_VORTEX(DM dm, SimpleBC bc, Physics phys,
206                                void *setup_ctx) {
207   PetscErrorCode ierr;
208   PetscFunctionBeginUser;
209 
210   // Define boundary conditions
211   bc->num_slip[2] = 2; bc->slips[2][0] = 1; bc->slips[2][1] = 2;
212 
213   // Set boundary conditions
214   DMLabel label;
215   ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr);
216   PetscInt comps[1] = {3};
217   ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, "Face Sets",
218                        bc->num_slip[2], bc->slips[2], 0, 1, comps,
219                        (void(*)(void))NULL, NULL, setup_ctx, NULL);
220   CHKERRQ(ierr);
221 
222   PetscFunctionReturn(0);
223 }
224 
225 PetscErrorCode PRINT_EULER_VORTEX(Physics phys, SetupContext setup_ctx,
226                                   AppCtx app_ctx) {
227   MPI_Comm       comm = PETSC_COMM_WORLD;
228   PetscErrorCode ierr;
229   PetscFunctionBeginUser;
230 
231   ierr = PetscPrintf(comm,
232                      "  Problem:\n"
233                      "    Problem Name                       : %s\n"
234                      "    Test Case                          : %s\n"
235                      "    Background Velocity                : %f,%f,%f\n"
236                      "    Stabilization                      : %s\n",
237                      app_ctx->problem_name, EulerTestTypes[phys->euler_test],
238                      phys->euler_ctx->mean_velocity[0],
239                      phys->euler_ctx->mean_velocity[1],
240                      phys->euler_ctx->mean_velocity[2],
241                      StabilizationTypes[phys->stab]); CHKERRQ(ierr);
242 
243   PetscFunctionReturn(0);
244 }
245