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_inflow = Advection_InOutFlow; 57 problem->apply_inflow_loc = Advection_InOutFlow_loc; 58 problem->bc = Exact_Advection; 59 problem->setup_ctx = SetupContext_ADVECTION; 60 problem->non_zero_time = PETSC_FALSE; 61 problem->print_info = PRINT_ADVECTION; 62 63 // ------------------------------------------------------ 64 // Create the libCEED context 65 // ------------------------------------------------------ 66 CeedScalar rc = 1000.; // m (Radius of bubble) 67 CeedScalar CtauS = 0.; // dimensionless 68 CeedScalar strong_form = 0.; // [0,1] 69 CeedScalar E_wind = 1.e6; // J 70 PetscReal wind[3] = {1., 0, 0}; // m/s 71 PetscReal domain_min[3], domain_max[3], domain_size[3]; 72 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 73 for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 74 75 76 // ------------------------------------------------------ 77 // Create the PETSc context 78 // ------------------------------------------------------ 79 PetscScalar meter = 1e-2; // 1 meter in scaled length units 80 PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units 81 PetscScalar second = 1e-2; // 1 second in scaled time units 82 PetscScalar Joule; 83 84 // ------------------------------------------------------ 85 // Command line Options 86 // ------------------------------------------------------ 87 ierr = PetscOptionsBegin(comm, NULL, "Options for ADVECTION problem", 88 NULL); CHKERRQ(ierr); 89 // -- Physics 90 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 91 NULL, rc, &rc, NULL); CHKERRQ(ierr); 92 PetscBool translation; 93 ierr = PetscOptionsEnum("-wind_type", "Wind type in Advection", 94 NULL, WindTypes, 95 (PetscEnum)(wind_type = WIND_ROTATION), 96 (PetscEnum *)&wind_type, &translation); CHKERRQ(ierr); 97 if (translation) user->phys->has_neumann = PETSC_TRUE; 98 PetscInt n = problem->dim; 99 PetscBool user_wind; 100 ierr = PetscOptionsRealArray("-wind_translation", "Constant wind vector", 101 NULL, wind, &n, &user_wind); CHKERRQ(ierr); 102 ierr = PetscOptionsScalar("-CtauS", 103 "Scale coefficient for tau (nondimensional)", 104 NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); 105 ierr = PetscOptionsScalar("-strong_form", 106 "Strong (1) or weak/integrated by parts (0) advection residual", 107 NULL, strong_form, &strong_form, NULL); CHKERRQ(ierr); 108 ierr = PetscOptionsScalar("-E_wind", "Total energy of inflow wind", 109 NULL, E_wind, &E_wind, NULL); CHKERRQ(ierr); 110 ierr = PetscOptionsEnum("-bubble_type", "Sphere (3D) or cylinder (2D)", 111 NULL, BubbleTypes, 112 (PetscEnum)(bubble_type = BUBBLE_SPHERE), 113 (PetscEnum *)&bubble_type, NULL); CHKERRQ(ierr); 114 ierr = PetscOptionsEnum("-bubble_continuity", "Smooth, back_sharp, or thick", 115 NULL, BubbleContinuityTypes, 116 (PetscEnum)(bubble_continuity_type = BUBBLE_CONTINUITY_SMOOTH), 117 (PetscEnum *)&bubble_continuity_type, NULL); CHKERRQ(ierr); 118 ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, 119 StabilizationTypes, (PetscEnum)(stab = STAB_NONE), 120 (PetscEnum *)&stab, NULL); CHKERRQ(ierr); 121 ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", 122 NULL, implicit=PETSC_FALSE, &implicit, NULL); 123 CHKERRQ(ierr); 124 125 // -- Units 126 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 127 NULL, meter, &meter, NULL); CHKERRQ(ierr); 128 meter = fabs(meter); 129 ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units", 130 NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr); 131 kilogram = fabs(kilogram); 132 ierr = PetscOptionsScalar("-units_second","1 second in scaled time units", 133 NULL, second, &second, NULL); CHKERRQ(ierr); 134 second = fabs(second); 135 136 // -- Warnings 137 if (wind_type == WIND_ROTATION && user_wind) { 138 ierr = PetscPrintf(comm, 139 "Warning! Use -wind_translation only with -wind_type translation\n"); 140 CHKERRQ(ierr); 141 } 142 if (wind_type == WIND_TRANSLATION 143 && bubble_type == BUBBLE_CYLINDER && wind[2] != 0.) { 144 wind[2] = 0; 145 ierr = PetscPrintf(comm, 146 "Warning! Background wind in the z direction should be zero (-wind_translation x,x,0) with -bubble_type cylinder\n"); 147 CHKERRQ(ierr); 148 } 149 if (stab == STAB_NONE && CtauS != 0) { 150 ierr = PetscPrintf(comm, 151 "Warning! Use -CtauS only with -stab su or -stab supg\n"); 152 CHKERRQ(ierr); 153 } 154 if (stab == STAB_SUPG && !implicit) { 155 ierr = PetscPrintf(comm, 156 "Warning! Use -stab supg only with -implicit\n"); 157 CHKERRQ(ierr); 158 } 159 160 ierr = PetscOptionsEnd(); CHKERRQ(ierr); 161 162 // ------------------------------------------------------ 163 // Set up the PETSc context 164 // ------------------------------------------------------ 165 // -- Define derived units 166 Joule = kilogram * PetscSqr(meter) / PetscSqr(second); 167 168 user->units->meter = meter; 169 user->units->kilogram = kilogram; 170 user->units->second = second; 171 user->units->Joule = Joule; 172 173 // ------------------------------------------------------ 174 // Set up the libCEED context 175 // ------------------------------------------------------ 176 // -- Scale variables to desired units 177 E_wind *= Joule; 178 rc = fabs(rc) * meter; 179 for (int i=0; i<3; i++) { 180 wind[i] *= (meter/second); 181 domain_size[i] *= meter; 182 } 183 problem->dm_scale = meter; 184 185 // -- Setup Context 186 setup_context->rc = rc; 187 setup_context->lx = domain_size[0]; 188 setup_context->ly = domain_size[1]; 189 setup_context->lz = domain_size[2]; 190 setup_context->wind[0] = wind[0]; 191 setup_context->wind[1] = wind[1]; 192 setup_context->wind[2] = wind[2]; 193 setup_context->wind_type = wind_type; 194 setup_context->bubble_type = bubble_type; 195 setup_context->bubble_continuity_type = bubble_continuity_type; 196 setup_context->time = 0; 197 198 // -- QFunction Context 199 user->phys->stab = stab; 200 user->phys->wind_type = wind_type; 201 user->phys->bubble_type = bubble_type; 202 user->phys->bubble_continuity_type = bubble_continuity_type; 203 // if passed correctly 204 user->phys->implicit = implicit; 205 user->phys->has_curr_time = has_curr_time; 206 user->phys->advection_ctx->CtauS = CtauS; 207 user->phys->advection_ctx->E_wind = E_wind; 208 user->phys->advection_ctx->implicit = implicit; 209 user->phys->advection_ctx->strong_form = strong_form; 210 user->phys->advection_ctx->stabilization = stab; 211 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode SetupContext_ADVECTION(Ceed ceed, CeedData ceed_data, 216 AppCtx app_ctx, SetupContext setup_ctx, Physics phys) { 217 PetscFunctionBeginUser; 218 CeedQFunctionContextCreate(ceed, &ceed_data->setup_context); 219 CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST, 220 CEED_USE_POINTER, sizeof(*setup_ctx), setup_ctx); 221 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->setup_context); 222 CeedQFunctionContextCreate(ceed, &ceed_data->advection_context); 223 CeedQFunctionContextSetData(ceed_data->advection_context, CEED_MEM_HOST, 224 CEED_USE_POINTER, 225 sizeof(*phys->advection_ctx), phys->advection_ctx); 226 if (ceed_data->qf_rhs_vol) 227 CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->advection_context); 228 if (ceed_data->qf_ifunction_vol) 229 CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, 230 ceed_data->advection_context); 231 if (ceed_data->qf_apply_inflow) 232 CeedQFunctionSetContext(ceed_data->qf_apply_inflow, 233 ceed_data->advection_context); 234 PetscFunctionReturn(0); 235 } 236 237 PetscErrorCode PRINT_ADVECTION(Physics phys, SetupContext setup_ctx, 238 AppCtx app_ctx) { 239 MPI_Comm comm = PETSC_COMM_WORLD; 240 PetscErrorCode ierr; 241 PetscFunctionBeginUser; 242 243 ierr = PetscPrintf(comm, 244 " Problem:\n" 245 " Problem Name : %s\n" 246 " Stabilization : %s\n" 247 " Bubble Type : %s (%dD)\n" 248 " Bubble Continuity : %s\n" 249 " Wind Type : %s\n", 250 app_ctx->problem_name, StabilizationTypes[phys->stab], 251 BubbleTypes[phys->bubble_type], 252 phys->bubble_type == BUBBLE_SPHERE ? 3 : 2, 253 BubbleContinuityTypes[phys->bubble_continuity_type], 254 WindTypes[phys->wind_type]); CHKERRQ(ierr); 255 256 if (phys->wind_type == WIND_TRANSLATION) { 257 ierr = PetscPrintf(comm, 258 " Background Wind : %f,%f,%f\n", 259 setup_ctx->wind[0], setup_ctx->wind[1], setup_ctx->wind[2]); CHKERRQ(ierr); 260 } 261 PetscFunctionReturn(0); 262 } 263