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