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