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