1 // Copyright (c) 2017-2022, 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 ADVECTION 10 11 #include "../navierstokes.h" 12 #include "../qfunctions/setupgeo.h" 13 #include "../qfunctions/advection.h" 14 15 PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, 16 void *ctx) { 17 WindType wind_type; 18 BubbleType bubble_type; 19 BubbleContinuityType bubble_continuity_type; 20 StabilizationType stab; 21 SetupContext setup_context; 22 User user = *(User *)ctx; 23 MPI_Comm comm = PETSC_COMM_WORLD; 24 PetscBool implicit; 25 PetscBool has_curr_time = PETSC_FALSE; 26 PetscInt ierr; 27 AdvectionContext advection_ctx; 28 CeedQFunctionContext advection_context; 29 30 PetscFunctionBeginUser; 31 ierr = PetscCalloc1(1, &setup_context); CHKERRQ(ierr); 32 ierr = PetscCalloc1(1, &advection_ctx); CHKERRQ(ierr); 33 34 // ------------------------------------------------------ 35 // SET UP ADVECTION 36 // ------------------------------------------------------ 37 problem->dim = 3; 38 problem->q_data_size_vol = 10; 39 problem->q_data_size_sur = 4; 40 problem->setup_vol.qfunction = Setup; 41 problem->setup_vol.qfunction_loc = Setup_loc; 42 problem->setup_sur.qfunction = SetupBoundary; 43 problem->setup_sur.qfunction_loc = SetupBoundary_loc; 44 problem->ics.qfunction = ICsAdvection; 45 problem->ics.qfunction_loc = ICsAdvection_loc; 46 problem->apply_vol_rhs.qfunction = Advection; 47 problem->apply_vol_rhs.qfunction_loc = Advection_loc; 48 problem->apply_vol_ifunction.qfunction = IFunction_Advection; 49 problem->apply_vol_ifunction.qfunction_loc = IFunction_Advection_loc; 50 problem->apply_inflow.qfunction = Advection_InOutFlow; 51 problem->apply_inflow.qfunction_loc = Advection_InOutFlow_loc; 52 problem->bc = Exact_Advection; 53 problem->bc_ctx = setup_context; 54 problem->non_zero_time = PETSC_FALSE; 55 problem->print_info = PRINT_ADVECTION; 56 57 // ------------------------------------------------------ 58 // Create the libCEED context 59 // ------------------------------------------------------ 60 CeedScalar rc = 1000.; // m (Radius of bubble) 61 CeedScalar CtauS = 0.; // dimensionless 62 CeedScalar strong_form = 0.; // [0,1] 63 CeedScalar E_wind = 1.e6; // J 64 PetscReal wind[3] = {1., 0, 0}; // m/s 65 PetscReal domain_min[3], domain_max[3], domain_size[3]; 66 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 67 for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 68 69 70 // ------------------------------------------------------ 71 // Create the PETSc context 72 // ------------------------------------------------------ 73 PetscScalar meter = 1e-2; // 1 meter in scaled length units 74 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 75 PetscScalar second = 1e-2; // 1 second in scaled time units 76 PetscScalar Joule; 77 78 // ------------------------------------------------------ 79 // Command line Options 80 // ------------------------------------------------------ 81 PetscOptionsBegin(comm, NULL, "Options for ADVECTION problem", NULL); 82 // -- Physics 83 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 84 NULL, rc, &rc, NULL); CHKERRQ(ierr); 85 PetscBool translation; 86 ierr = PetscOptionsEnum("-wind_type", "Wind type in Advection", 87 NULL, WindTypes, 88 (PetscEnum)(wind_type = WIND_ROTATION), 89 (PetscEnum *)&wind_type, &translation); CHKERRQ(ierr); 90 if (translation) user->phys->has_neumann = PETSC_TRUE; 91 PetscInt n = problem->dim; 92 PetscBool user_wind; 93 ierr = PetscOptionsRealArray("-wind_translation", "Constant wind vector", 94 NULL, wind, &n, &user_wind); CHKERRQ(ierr); 95 ierr = PetscOptionsScalar("-CtauS", 96 "Scale coefficient for tau (nondimensional)", 97 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 98 ierr = PetscOptionsScalar("-strong_form", 99 "Strong (1) or weak/integrated by parts (0) advection residual", 100 NULL, strong_form, &strong_form, NULL); CHKERRQ(ierr); 101 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 102 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 103 ierr = PetscOptionsEnum("-bubble_type", "Sphere (3D) or cylinder (2D)", 104 NULL, BubbleTypes, 105 (PetscEnum)(bubble_type = BUBBLE_SPHERE), 106 (PetscEnum *)&bubble_type, NULL); CHKERRQ(ierr); 107 ierr = PetscOptionsEnum("-bubble_continuity", "Smooth, back_sharp, or thick", 108 NULL, BubbleContinuityTypes, 109 (PetscEnum)(bubble_continuity_type = BUBBLE_CONTINUITY_SMOOTH), 110 (PetscEnum *)&bubble_continuity_type, NULL); CHKERRQ(ierr); 111 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 112 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 113 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 114 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 115 NULL, implicit=PETSC_FALSE, &implicit, NULL); 116 CHKERRQ(ierr); 117 118 // -- Units 119 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 120 NULL, meter, &meter, NULL); CHKERRQ(ierr); 121 meter = fabs(meter); 122 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 123 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 124 kilogram = fabs(kilogram); 125 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 126 NULL, second, &second, NULL); CHKERRQ(ierr); 127 second = fabs(second); 128 129 // -- Warnings 130 if (wind_type == WIND_ROTATION && user_wind) { 131 ierr = PetscPrintf(comm, 132 "Warning! Use -wind_translation only with -wind_type translation\n"); 133 CHKERRQ(ierr); 134 } 135 if (wind_type == WIND_TRANSLATION 136 && bubble_type == BUBBLE_CYLINDER && wind[2] != 0.) { 137 wind[2] = 0; 138 ierr = PetscPrintf(comm, 139 "Warning! Background wind in the z direction should be zero (-wind_translation x,x,0) with -bubble_type cylinder\n"); 140 CHKERRQ(ierr); 141 } 142 if (stab == STAB_NONE && CtauS != 0) { 143 ierr = PetscPrintf(comm, 144 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 145 CHKERRQ(ierr); 146 } 147 if (stab == STAB_SUPG && !implicit) { 148 ierr = PetscPrintf(comm, 149 "Warning! Use -stab supg only with -implicit\n"); 150 CHKERRQ(ierr); 151 } 152 153 PetscOptionsEnd(); 154 155 // ------------------------------------------------------ 156 // Set up the PETSc context 157 // ------------------------------------------------------ 158 // -- Define derived units 159 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 160 161 user->units->meter = meter; 162 user->units->kilogram = kilogram; 163 user->units->second = second; 164 user->units->Joule = Joule; 165 166 // ------------------------------------------------------ 167 // Set up the libCEED context 168 // ------------------------------------------------------ 169 // -- Scale variables to desired units 170 E_wind *= Joule; 171 rc = fabs(rc) * meter; 172 for (int i=0; i<3; i++) { 173 wind[i] *= (meter/second); 174 domain_size[i] *= meter; 175 } 176 problem->dm_scale = meter; 177 178 // -- Setup Context 179 setup_context->rc = rc; 180 setup_context->lx = domain_size[0]; 181 setup_context->ly = domain_size[1]; 182 setup_context->lz = domain_size[2]; 183 setup_context->wind[0] = wind[0]; 184 setup_context->wind[1] = wind[1]; 185 setup_context->wind[2] = wind[2]; 186 setup_context->wind_type = wind_type; 187 setup_context->bubble_type = bubble_type; 188 setup_context->bubble_continuity_type = bubble_continuity_type; 189 setup_context->time = 0; 190 191 // -- QFunction Context 192 user->phys->stab = stab; 193 user->phys->wind_type = wind_type; 194 user->phys->bubble_type = bubble_type; 195 user->phys->bubble_continuity_type = bubble_continuity_type; 196 // if passed correctly 197 user->phys->implicit = implicit; 198 user->phys->has_curr_time = has_curr_time; 199 advection_ctx->CtauS = CtauS; 200 advection_ctx->E_wind = E_wind; 201 advection_ctx->implicit = implicit; 202 advection_ctx->strong_form = strong_form; 203 advection_ctx->stabilization = stab; 204 205 CeedQFunctionContextCreate(user->ceed, &problem->ics.qfunction_context); 206 CeedQFunctionContextSetData(problem->ics.qfunction_context, CEED_MEM_HOST, 207 CEED_USE_POINTER, sizeof(*setup_context), setup_context); 208 209 CeedQFunctionContextCreate(user->ceed, &advection_context); 210 CeedQFunctionContextSetData(advection_context, CEED_MEM_HOST, 211 CEED_USE_POINTER, 212 sizeof(*advection_ctx), advection_ctx); 213 CeedQFunctionContextSetDataDestroy(advection_context, CEED_MEM_HOST, 214 FreeContextPetsc); 215 problem->apply_vol_rhs.qfunction_context = advection_context; 216 CeedQFunctionContextReferenceCopy(advection_context, 217 &problem->apply_vol_ifunction.qfunction_context); 218 CeedQFunctionContextReferenceCopy(advection_context, 219 &problem->apply_inflow.qfunction_context); 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode PRINT_ADVECTION(ProblemData *problem, AppCtx app_ctx) { 224 MPI_Comm comm = PETSC_COMM_WORLD; 225 PetscErrorCode ierr; 226 SetupContext setup_ctx; 227 AdvectionContext advection_ctx; 228 229 PetscFunctionBeginUser; 230 CeedQFunctionContextGetData(problem->ics.qfunction_context, 231 CEED_MEM_HOST, &setup_ctx); 232 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 233 CEED_MEM_HOST, &advection_ctx); 234 ierr = PetscPrintf(comm, 235 " Problem:\n" 236 " Problem Name : %s\n" 237 " Stabilization : %s\n" 238 " Bubble Type : %s (%dD)\n" 239 " Bubble Continuity : %s\n" 240 " Wind Type : %s\n", 241 app_ctx->problem_name, StabilizationTypes[advection_ctx->stabilization], 242 BubbleTypes[setup_ctx->bubble_type], 243 setup_ctx->bubble_type == BUBBLE_SPHERE ? 3 : 2, 244 BubbleContinuityTypes[setup_ctx->bubble_continuity_type], 245 WindTypes[setup_ctx->wind_type]); CHKERRQ(ierr); 246 247 if (setup_ctx->wind_type == WIND_TRANSLATION) { 248 ierr = PetscPrintf(comm, 249 " Background Wind : %f,%f,%f\n", 250 setup_ctx->wind[0], setup_ctx->wind[1], setup_ctx->wind[2]); CHKERRQ(ierr); 251 } 252 CeedQFunctionContextRestoreData(problem->ics.qfunction_context, 253 &setup_ctx); 254 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 255 &advection_ctx); 256 PetscFunctionReturn(0); 257 } 258