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 Channel flow 10 11 #include "../navierstokes.h" 12 #include "../qfunctions/channel.h" 13 14 PetscErrorCode NS_CHANNEL(ProblemData *problem, DM dm, void *ctx) { 15 16 PetscInt ierr; 17 User user = *(User *)ctx; 18 MPI_Comm comm = PETSC_COMM_WORLD; 19 ChannelContext channel_ctx; 20 NewtonianIdealGasContext newtonian_ig_ctx; 21 CeedQFunctionContext channel_context; 22 23 PetscFunctionBeginUser; 24 ierr = NS_NEWTONIAN_IG(problem, dm, ctx); CHKERRQ(ierr); 25 ierr = PetscCalloc1(1, &channel_ctx); CHKERRQ(ierr); 26 27 // ------------------------------------------------------ 28 // SET UP Channel 29 // ------------------------------------------------------ 30 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 31 problem->ics.qfunction = ICsChannel; 32 problem->ics.qfunction_loc = ICsChannel_loc; 33 if (user->phys->state_var == STATEVAR_CONSERVATIVE) { 34 problem->apply_inflow.qfunction = Channel_Inflow; 35 problem->apply_inflow.qfunction_loc = Channel_Inflow_loc; 36 problem->apply_outflow.qfunction = Channel_Outflow; 37 problem->apply_outflow.qfunction_loc = Channel_Outflow_loc; 38 } 39 40 // -- Command Line Options 41 CeedScalar umax = 10.; // m/s 42 CeedScalar theta0 = 300.; // K 43 CeedScalar P0 = 1.e5; // Pa 44 PetscReal body_force_scale = 1.; 45 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL); 46 ierr = PetscOptionsScalar("-umax", "Centerline velocity of the Channel", 47 NULL, umax, &umax, NULL); CHKERRQ(ierr); 48 ierr = PetscOptionsScalar("-theta0", "Wall temperature", 49 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 50 ierr = PetscOptionsScalar("-P0", "Pressure at outflow", 51 NULL, P0, &P0, NULL); CHKERRQ(ierr); 52 ierr = PetscOptionsReal("-body_force_scale", "Multiplier for body force", 53 NULL, body_force_scale=1, &body_force_scale, NULL); CHKERRQ(ierr); 54 PetscOptionsEnd(); 55 56 PetscScalar meter = user->units->meter; 57 PetscScalar second = user->units->second; 58 PetscScalar Kelvin = user->units->Kelvin; 59 PetscScalar Pascal = user->units->Pascal; 60 61 theta0 *= Kelvin; 62 P0 *= Pascal; 63 umax *= meter / second; 64 65 //-- Setup Problem information 66 CeedScalar H, center; 67 { 68 PetscReal domain_min[3], domain_max[3], domain_size[3]; 69 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 70 for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 71 72 H = 0.5*domain_size[1]*meter; 73 center = H + domain_min[1]*meter; 74 } 75 76 // Some properties depend on parameters from NewtonianIdealGas 77 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 78 CEED_MEM_HOST, &newtonian_ig_ctx); 79 80 channel_ctx->center = center; 81 channel_ctx->H = H; 82 channel_ctx->theta0 = theta0; 83 channel_ctx->P0 = P0; 84 channel_ctx->umax = umax; 85 channel_ctx->implicit = user->phys->implicit; 86 channel_ctx->B = body_force_scale * 2 * umax*newtonian_ig_ctx->mu / (H*H); 87 88 { 89 // Calculate Body force 90 CeedScalar cv = newtonian_ig_ctx->cv, 91 cp = newtonian_ig_ctx->cp; 92 CeedScalar Rd = cp - cv; 93 CeedScalar rho = P0 / (Rd*theta0); 94 CeedScalar g[] = {channel_ctx->B / rho, 0., 0.}; 95 ierr = PetscArraycpy(newtonian_ig_ctx->g, g, 3); CHKERRQ(ierr); 96 } 97 channel_ctx->newtonian_ctx = *newtonian_ig_ctx; 98 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 99 &newtonian_ig_ctx); 100 101 CeedQFunctionContextCreate(user->ceed, &channel_context); 102 CeedQFunctionContextSetData(channel_context, CEED_MEM_HOST, 103 CEED_USE_POINTER, 104 sizeof(*channel_ctx), channel_ctx); 105 CeedQFunctionContextSetDataDestroy(channel_context, CEED_MEM_HOST, 106 FreeContextPetsc); 107 108 problem->ics.qfunction_context = channel_context; 109 CeedQFunctionContextReferenceCopy(channel_context, 110 &problem->apply_inflow.qfunction_context); 111 CeedQFunctionContextReferenceCopy(channel_context, 112 &problem->apply_outflow.qfunction_context); 113 PetscFunctionReturn(0); 114 } 115