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