1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Utility functions for setting up DENSITY_CURRENT 6 7 #include "../qfunctions/densitycurrent.h" 8 9 #include <ceed.h> 10 #include <petscdm.h> 11 12 #include <navierstokes.h> 13 14 PetscErrorCode NS_DENSITY_CURRENT(ProblemData problem, DM dm, void *ctx) { 15 Honee honee = *(Honee *)ctx; 16 MPI_Comm comm = honee->comm; 17 Ceed ceed = honee->ceed; 18 DensityCurrentContext density_current_ctx; 19 CeedQFunctionContext density_current_qfctx; 20 NewtonianIdealGasContext newtonian_ig_ctx; 21 PetscInt dim; 22 23 PetscFunctionBeginUser; 24 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx)); 25 PetscCall(PetscNew(&density_current_ctx)); 26 27 // ------------------------------------------------------ 28 // Create the QFunction context 29 // ------------------------------------------------------ 30 CeedScalar theta0 = 300.; // K 31 CeedScalar thetaC = -15.; // K 32 CeedScalar P0 = 1.e5; // Pa 33 CeedScalar N = 0.01; // 1/s 34 CeedScalar rc = 1000.; // m (Radius of bubble) 35 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 36 PetscReal domain_min[3], domain_max[3], domain_size[3]; 37 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 38 PetscCall(DMGetDimension(dm, &dim)); 39 for (PetscInt i = 0; i < dim; i++) domain_size[i] = domain_max[i] - domain_min[i]; 40 41 // ------------------------------------------------------ 42 // Command line Options 43 // ------------------------------------------------------ 44 PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem", NULL); 45 PetscCall(PetscOptionsScalar("-theta0", "Reference potential temperature", NULL, theta0, &theta0, NULL)); 46 PetscCall(PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", NULL, thetaC, &thetaC, NULL)); 47 PetscCall(PetscOptionsScalar("-P0", "Atmospheric pressure", NULL, P0, &P0, NULL)); 48 PetscCall(PetscOptionsScalar("-N", "Brunt-Vaisala frequency", NULL, N, &N, NULL)); 49 PetscCall(PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", NULL, rc, &rc, NULL)); 50 for (PetscInt i = 0; i < dim; i++) center[i] = .5 * domain_size[i]; 51 PetscInt n = dim; 52 PetscCall(PetscOptionsRealArray("-center", "Location of bubble center", NULL, center, &n, NULL)); 53 n = dim; 54 PetscCall(PetscOptionsRealArray("-dc_axis", 55 "Axis of density current cylindrical anomaly, " 56 "or {0,0,0} for spherically symmetric", 57 NULL, dc_axis, &n, NULL)); 58 { 59 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 60 if (norm > 0) { 61 for (PetscInt i = 0; i < dim; i++) dc_axis[i] /= norm; 62 } 63 } 64 65 PetscOptionsEnd(); 66 67 Units units = honee->units; 68 69 rc = fabs(rc) * units->meter; 70 theta0 *= units->Kelvin; 71 thetaC *= units->Kelvin; 72 P0 *= units->Pascal; 73 N *= (1. / units->second); 74 for (PetscInt i = 0; i < dim; i++) center[i] *= units->meter; 75 76 density_current_ctx->theta0 = theta0; 77 density_current_ctx->thetaC = thetaC; 78 density_current_ctx->P0 = P0; 79 density_current_ctx->N = N; 80 density_current_ctx->rc = rc; 81 PetscCall(PetscArraycpy(density_current_ctx->center, center, 3)); 82 PetscCall(PetscArraycpy(density_current_ctx->dc_axis, dc_axis, 3)); 83 84 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 85 density_current_ctx->newtonian_ctx = *newtonian_ig_ctx; 86 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 87 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &density_current_qfctx)); 88 PetscCallCeed( 89 ceed, CeedQFunctionContextSetData(density_current_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*density_current_ctx), density_current_ctx)); 90 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(density_current_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 91 92 // ------------------------------------------------------ 93 // SET UP DENSITY_CURRENT 94 // ------------------------------------------------------ 95 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 96 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsDC, .qf_loc = ICsDC_loc, .qfctx = density_current_qfctx}; 97 PetscFunctionReturn(PETSC_SUCCESS); 98 } 99