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