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 DENSITY_CURRENT 10 11 #include "../qfunctions/densitycurrent.h" 12 13 #include <ceed.h> 14 #include <petscdm.h> 15 16 #include "../navierstokes.h" 17 18 PetscErrorCode NS_DENSITY_CURRENT(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 DensityCurrentContext dc_ctx; 23 CeedQFunctionContext density_current_context; 24 NewtonianIdealGasContext newtonian_ig_ctx; 25 26 PetscFunctionBeginUser; 27 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx, bc)); 28 PetscCall(PetscCalloc1(1, &dc_ctx)); 29 // ------------------------------------------------------ 30 // SET UP DENSITY_CURRENT 31 // ------------------------------------------------------ 32 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context)); 33 problem->ics.qfunction = ICsDC; 34 problem->ics.qfunction_loc = ICsDC_loc; 35 36 // ------------------------------------------------------ 37 // Create the libCEED context 38 // ------------------------------------------------------ 39 CeedScalar theta0 = 300.; // K 40 CeedScalar thetaC = -15.; // K 41 CeedScalar P0 = 1.e5; // Pa 42 CeedScalar N = 0.01; // 1/s 43 CeedScalar rc = 1000.; // m (Radius of bubble) 44 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 45 PetscReal domain_min[3], domain_max[3], domain_size[3]; 46 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 47 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 48 49 // ------------------------------------------------------ 50 // Command line Options 51 // ------------------------------------------------------ 52 PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem", NULL); 53 PetscCall(PetscOptionsScalar("-theta0", "Reference potential temperature", NULL, theta0, &theta0, NULL)); 54 PetscCall(PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", NULL, thetaC, &thetaC, NULL)); 55 PetscCall(PetscOptionsScalar("-P0", "Atmospheric pressure", NULL, P0, &P0, NULL)); 56 PetscCall(PetscOptionsScalar("-N", "Brunt-Vaisala frequency", NULL, N, &N, NULL)); 57 PetscCall(PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", NULL, rc, &rc, NULL)); 58 for (PetscInt i = 0; i < 3; i++) center[i] = .5 * domain_size[i]; 59 PetscInt n = problem->dim; 60 PetscCall(PetscOptionsRealArray("-center", "Location of bubble center", NULL, center, &n, NULL)); 61 n = problem->dim; 62 PetscCall(PetscOptionsRealArray("-dc_axis", 63 "Axis of density current cylindrical anomaly, " 64 "or {0,0,0} for spherically symmetric", 65 NULL, dc_axis, &n, NULL)); 66 { 67 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 68 if (norm > 0) { 69 for (PetscInt i = 0; i < 3; i++) dc_axis[i] /= norm; 70 } 71 } 72 73 PetscOptionsEnd(); 74 75 PetscScalar meter = user->units->meter; 76 PetscScalar second = user->units->second; 77 PetscScalar Kelvin = user->units->Kelvin; 78 PetscScalar Pascal = user->units->Pascal; 79 rc = fabs(rc) * meter; 80 theta0 *= Kelvin; 81 thetaC *= Kelvin; 82 P0 *= Pascal; 83 N *= (1. / second); 84 for (PetscInt i = 0; i < 3; i++) center[i] *= meter; 85 86 dc_ctx->theta0 = theta0; 87 dc_ctx->thetaC = thetaC; 88 dc_ctx->P0 = P0; 89 dc_ctx->N = N; 90 dc_ctx->rc = rc; 91 dc_ctx->center[0] = center[0]; 92 dc_ctx->center[1] = center[1]; 93 dc_ctx->center[2] = center[2]; 94 dc_ctx->dc_axis[0] = dc_axis[0]; 95 dc_ctx->dc_axis[1] = dc_axis[1]; 96 dc_ctx->dc_axis[2] = dc_axis[2]; 97 98 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx)); 99 dc_ctx->newtonian_ctx = *newtonian_ig_ctx; 100 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx)); 101 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &density_current_context)); 102 PetscCallCeed(ceed, CeedQFunctionContextSetData(density_current_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*dc_ctx), dc_ctx)); 103 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(density_current_context, CEED_MEM_HOST, FreeContextPetsc)); 104 problem->ics.qfunction_context = density_current_context; 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107