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 Blasius Boundary Layer 10 11 #include "../navierstokes.h" 12 #include "../qfunctions/blasius.h" 13 #include "stg_shur14.h" 14 15 /* \brief Modify the domain and mesh for blasius 16 * 17 * Modifies mesh such that `N` elements are within `refine_height` with a 18 * geometric growth ratio of `growth`. Excess elements are then distributed 19 * linearly in logspace to the top surface. 20 * 21 * The top surface is also angled downwards, so that it may be used as an 22 * outflow. It's angle is controlled by `top_angle` (in units of degrees). 23 */ 24 static PetscErrorCode ModifyMesh(DM dm, PetscInt dim, PetscReal growth, 25 PetscInt N, PetscReal refine_height, 26 PetscReal top_angle) { 27 28 PetscInt ierr, narr, ncoords; 29 PetscReal domain_min[3], domain_max[3], domain_size[3]; 30 PetscScalar *arr_coords; 31 Vec vec_coords; 32 PetscFunctionBeginUser; 33 34 PetscReal angle_coeff = tan(top_angle*(M_PI/180)); 35 36 // Get domain boundary information 37 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 38 for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 39 40 // Get coords array from DM 41 ierr = DMGetCoordinatesLocal(dm, &vec_coords); CHKERRQ(ierr); 42 ierr = VecGetLocalSize(vec_coords, &narr); CHKERRQ(ierr); 43 ierr = VecGetArray(vec_coords, &arr_coords); CHKERRQ(ierr); 44 45 PetscScalar (*coords)[dim] = (PetscScalar(*)[dim]) arr_coords; 46 ncoords = narr/dim; 47 48 // Get mesh information 49 PetscInt nmax = 3, faces[3]; 50 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, 51 NULL); CHKERRQ(ierr); 52 53 // Calculate the first element height 54 PetscReal dybox = domain_size[1]/faces[1]; 55 PetscReal dy1 = refine_height*(growth-1)/(pow(growth, N)-1); 56 57 // Calculate log of sizing outside BL 58 PetscReal logdy = (log(domain_max[1]) - log(refine_height)) / (faces[1] - N); 59 60 for(PetscInt i=0; i<ncoords; i++) { 61 PetscInt y_box_index = round(coords[i][1]/dybox); 62 if(y_box_index <= N) { 63 coords[i][1] = (1 - (coords[i][0]/domain_max[0])*angle_coeff) * 64 dy1*(pow(growth, coords[i][1]/dybox)-1)/(growth-1); 65 } else { 66 PetscInt j = y_box_index - N; 67 coords[i][1] = (1 - (coords[i][0]/domain_max[0])*angle_coeff) * 68 exp(log(refine_height) + logdy*j); 69 } 70 } 71 72 ierr = VecRestoreArray(vec_coords, &arr_coords); CHKERRQ(ierr); 73 ierr = DMSetCoordinatesLocal(dm, vec_coords); CHKERRQ(ierr); 74 75 PetscFunctionReturn(0); 76 } 77 78 PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx) { 79 80 PetscInt ierr; 81 User user = *(User *)ctx; 82 MPI_Comm comm = PETSC_COMM_WORLD; 83 PetscBool use_stg = PETSC_FALSE; 84 BlasiusContext blasius_ctx; 85 NewtonianIdealGasContext newtonian_ig_ctx; 86 CeedQFunctionContext blasius_context; 87 88 PetscFunctionBeginUser; 89 ierr = NS_NEWTONIAN_IG(problem, dm, ctx); CHKERRQ(ierr); 90 ierr = PetscCalloc1(1, &blasius_ctx); CHKERRQ(ierr); 91 92 // ------------------------------------------------------ 93 // SET UP Blasius 94 // ------------------------------------------------------ 95 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 96 problem->ics.qfunction = ICsBlasius; 97 problem->ics.qfunction_loc = ICsBlasius_loc; 98 problem->apply_outflow.qfunction = Blasius_Outflow; 99 problem->apply_outflow.qfunction_loc = Blasius_Outflow_loc; 100 problem->apply_inflow.qfunction = Blasius_Inflow; 101 problem->apply_inflow.qfunction_loc = Blasius_Inflow_loc; 102 103 CeedScalar Uinf = 40; // m/s 104 CeedScalar delta0 = 4.2e-4; // m 105 CeedScalar theta0 = 288.; // K 106 CeedScalar P0 = 1.01e5; // Pa 107 PetscBool weakT = PETSC_FALSE; // weak density or temperature 108 PetscReal mesh_refine_height = 5.9e-4; // m 109 PetscReal mesh_growth = 1.08; // [-] 110 PetscInt mesh_Ndelta = 45; // [-] 111 PetscReal mesh_top_angle = 5; // degrees 112 113 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL); 114 ierr = PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", 115 NULL, weakT, &weakT, NULL); CHKERRQ(ierr); 116 ierr = PetscOptionsScalar("-Uinf", "Velocity at boundary layer edge", 117 NULL, Uinf, &Uinf, NULL); CHKERRQ(ierr); 118 ierr = PetscOptionsScalar("-delta0", "Boundary layer height at inflow", 119 NULL, delta0, &delta0, NULL); CHKERRQ(ierr); 120 ierr = PetscOptionsScalar("-theta0", "Wall temperature", 121 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 122 ierr = PetscOptionsScalar("-P0", "Pressure at outflow", 123 NULL, P0, &P0, NULL); CHKERRQ(ierr); 124 ierr = PetscOptionsBoundedInt("-platemesh_Ndelta", 125 "Velocity at boundary layer edge", 126 NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1); CHKERRQ(ierr); 127 ierr = PetscOptionsScalar("-platemesh_refine_height", 128 "Height of boundary layer mesh refinement", 129 NULL, mesh_refine_height, &mesh_refine_height, NULL); CHKERRQ(ierr); 130 ierr = PetscOptionsScalar("-platemesh_growth", 131 "Geometric growth rate of boundary layer mesh", 132 NULL, mesh_growth, &mesh_growth, NULL); CHKERRQ(ierr); 133 ierr = PetscOptionsScalar("-platemesh_top_angle", 134 "Geometric top_angle rate of boundary layer mesh", 135 NULL, mesh_top_angle, &mesh_top_angle, NULL); CHKERRQ(ierr); 136 ierr = PetscOptionsBool("-stg_use", "Use STG inflow boundary condition", 137 NULL, use_stg, &use_stg, NULL); CHKERRQ(ierr); 138 PetscOptionsEnd(); 139 140 PetscScalar meter = user->units->meter; 141 PetscScalar second = user->units->second; 142 PetscScalar Kelvin = user->units->Kelvin; 143 PetscScalar Pascal = user->units->Pascal; 144 145 theta0 *= Kelvin; 146 P0 *= Pascal; 147 Uinf *= meter / second; 148 delta0 *= meter; 149 150 ierr = ModifyMesh(dm, problem->dim, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle); 151 CHKERRQ(ierr); 152 153 // Some properties depend on parameters from NewtonianIdealGas 154 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 155 CEED_MEM_HOST, &newtonian_ig_ctx); 156 157 blasius_ctx->weakT = weakT; 158 blasius_ctx->Uinf = Uinf; 159 blasius_ctx->delta0 = delta0; 160 blasius_ctx->theta0 = theta0; 161 blasius_ctx->P0 = P0; 162 blasius_ctx->implicit = user->phys->implicit; 163 blasius_ctx->newtonian_ctx = *newtonian_ig_ctx; 164 165 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 166 &newtonian_ig_ctx); 167 168 CeedQFunctionContextCreate(user->ceed, &blasius_context); 169 CeedQFunctionContextSetData(blasius_context, CEED_MEM_HOST, 170 CEED_USE_POINTER, 171 sizeof(*blasius_ctx), blasius_ctx); 172 CeedQFunctionContextSetDataDestroy(blasius_context, CEED_MEM_HOST, 173 FreeContextPetsc); 174 175 problem->ics.qfunction_context = blasius_context; 176 CeedQFunctionContextReferenceCopy(blasius_context, 177 &problem->apply_inflow.qfunction_context); 178 CeedQFunctionContextReferenceCopy(blasius_context, 179 &problem->apply_outflow.qfunction_context); 180 if (use_stg) { 181 ierr = SetupSTG(comm, dm, problem, user, weakT, theta0, P0); CHKERRQ(ierr); 182 } 183 PetscFunctionReturn(0); 184 } 185