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 #ifndef channel_context_struct 15 #define channel_context_struct 16 typedef struct ChannelContext_ *ChannelContext; 17 struct ChannelContext_ { 18 bool implicit; // !< Using implicit timesteping or not 19 CeedScalar theta0; // !< Reference temperature 20 CeedScalar P0; // !< Reference Pressure 21 CeedScalar umax; // !< Centerline velocity 22 CeedScalar center; // !< Y Coordinate for center of channel 23 CeedScalar H; // !< Channel half-height 24 CeedScalar B; // !< Body-force driving the flow 25 struct NewtonianIdealGasContext_ newtonian_ctx; 26 }; 27 #endif 28 29 PetscErrorCode NS_CHANNEL(ProblemData *problem, DM dm, void *setup_ctx, 30 void *ctx) { 31 32 PetscInt ierr; 33 ierr = NS_NEWTONIAN_IG(problem, dm, setup_ctx, ctx); CHKERRQ(ierr); 34 User user = *(User *)ctx; 35 MPI_Comm comm = PETSC_COMM_WORLD; 36 PetscFunctionBeginUser; 37 ierr = PetscCalloc1(1, &user->phys->channel_ctx); CHKERRQ(ierr); 38 39 // ------------------------------------------------------ 40 // SET UP Channel 41 // ------------------------------------------------------ 42 problem->ics = ICsChannel; 43 problem->ics_loc = ICsChannel_loc; 44 problem->apply_inflow = Channel_Inflow; 45 problem->apply_inflow_loc = Channel_Inflow_loc; 46 problem->apply_outflow = Channel_Outflow; 47 problem->apply_outflow_loc = Channel_Outflow_loc; 48 problem->setup_ctx = SetupContext_CHANNEL; 49 50 // -- Command Line Options 51 CeedScalar umax = 10.; // m/s 52 CeedScalar mu = .01; // Pa s, dynamic viscosity 53 //TODO ^^ make optional/respect explicit user set 54 CeedScalar theta0 = 300.; // K 55 CeedScalar P0 = 1.e5; // Pa 56 PetscOptionsBegin(comm, NULL, "Options for CHANNEL problem", NULL); 57 ierr = PetscOptionsScalar("-umax", "Centerline velocity of the Channel", 58 NULL, umax, &umax, NULL); CHKERRQ(ierr); 59 ierr = PetscOptionsScalar("-theta0", "Wall temperature", 60 NULL, theta0, &theta0, NULL); CHKERRQ(ierr); 61 ierr = PetscOptionsScalar("-P0", "Pressure at outflow", 62 NULL, P0, &P0, NULL); CHKERRQ(ierr); 63 PetscOptionsEnd(); 64 65 PetscScalar meter = user->units->meter; 66 PetscScalar second = user->units->second; 67 PetscScalar Kelvin = user->units->Kelvin; 68 PetscScalar Pascal = user->units->Pascal; 69 70 mu *= Pascal * second; 71 theta0 *= Kelvin; 72 P0 *= Pascal; 73 umax *= meter / second; 74 75 //-- Setup Problem information 76 CeedScalar H, center; 77 { 78 PetscReal domain_min[3], domain_max[3], domain_size[3]; 79 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 80 for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 81 82 H = 0.5*domain_size[1]*meter; 83 center = H + domain_min[1]*meter; 84 } 85 86 user->phys->channel_ctx->center = center; 87 user->phys->channel_ctx->H = H; 88 user->phys->channel_ctx->theta0 = theta0; 89 user->phys->channel_ctx->P0 = P0; 90 user->phys->channel_ctx->umax = umax; 91 user->phys->channel_ctx->implicit = user->phys->implicit; 92 user->phys->channel_ctx->B = -2*umax*mu/H; 93 94 { 95 // Calculate Body force 96 CeedScalar cv = user->phys->newtonian_ig_ctx->cv, 97 cp = user->phys->newtonian_ig_ctx->cp; 98 CeedScalar Rd = cp - cv; 99 CeedScalar rho = P0 / (Rd*theta0); 100 CeedScalar g[] = {user->phys->channel_ctx->B / rho, 0., 0.}; 101 ierr = PetscArraycpy(user->phys->newtonian_ig_ctx->g, g, 3); CHKERRQ(ierr); 102 } 103 user->phys->newtonian_ig_ctx->mu = mu; 104 user->phys->channel_ctx->newtonian_ctx = *user->phys->newtonian_ig_ctx; 105 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode SetupContext_CHANNEL(Ceed ceed, CeedData ceed_data, 110 AppCtx app_ctx, SetupContext setup_ctx, Physics phys) { 111 PetscFunctionBeginUser; 112 PetscInt ierr; 113 ierr = SetupContext_NEWTONIAN_IG(ceed, ceed_data, app_ctx, setup_ctx, phys); 114 CHKERRQ(ierr); 115 CeedQFunctionContextCreate(ceed, &ceed_data->channel_context); 116 CeedQFunctionContextSetData(ceed_data->channel_context, CEED_MEM_HOST, 117 CEED_USE_POINTER, 118 sizeof(*phys->channel_ctx), phys->channel_ctx); 119 phys->has_neumann = PETSC_TRUE; 120 if (ceed_data->qf_ics) 121 CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->channel_context); 122 if (ceed_data->qf_apply_inflow) 123 CeedQFunctionSetContext(ceed_data->qf_apply_inflow, ceed_data->channel_context); 124 if (ceed_data->qf_apply_outflow) 125 CeedQFunctionSetContext(ceed_data->qf_apply_outflow, 126 ceed_data->channel_context); 127 PetscFunctionReturn(0); 128 } 129