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 // SET UP DENSITY_CURRENT 28 // ------------------------------------------------------ 29 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 30 problem->ics.qf_func_ptr = ICsDC; 31 problem->ics.qf_loc = ICsDC_loc; 32 33 // ------------------------------------------------------ 34 // Create the QFunction context 35 // ------------------------------------------------------ 36 CeedScalar theta0 = 300.; // K 37 CeedScalar thetaC = -15.; // K 38 CeedScalar P0 = 1.e5; // Pa 39 CeedScalar N = 0.01; // 1/s 40 CeedScalar rc = 1000.; // m (Radius of bubble) 41 PetscReal center[3], dc_axis[3] = {0, 0, 0}; 42 PetscReal domain_min[3], domain_max[3], domain_size[3]; 43 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 44 PetscCall(DMGetDimension(dm, &dim)); 45 for (PetscInt i = 0; i < dim; i++) domain_size[i] = domain_max[i] - domain_min[i]; 46 47 // ------------------------------------------------------ 48 // Command line Options 49 // ------------------------------------------------------ 50 PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem", NULL); 51 PetscCall(PetscOptionsScalar("-theta0", "Reference potential temperature", NULL, theta0, &theta0, NULL)); 52 PetscCall(PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", NULL, thetaC, &thetaC, NULL)); 53 PetscCall(PetscOptionsScalar("-P0", "Atmospheric pressure", NULL, P0, &P0, NULL)); 54 PetscCall(PetscOptionsScalar("-N", "Brunt-Vaisala frequency", NULL, N, &N, NULL)); 55 PetscCall(PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", NULL, rc, &rc, NULL)); 56 for (PetscInt i = 0; i < dim; i++) center[i] = .5 * domain_size[i]; 57 PetscInt n = dim; 58 PetscCall(PetscOptionsRealArray("-center", "Location of bubble center", NULL, center, &n, NULL)); 59 n = dim; 60 PetscCall(PetscOptionsRealArray("-dc_axis", 61 "Axis of density current cylindrical anomaly, " 62 "or {0,0,0} for spherically symmetric", 63 NULL, dc_axis, &n, NULL)); 64 { 65 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2])); 66 if (norm > 0) { 67 for (PetscInt i = 0; i < dim; i++) dc_axis[i] /= norm; 68 } 69 } 70 71 PetscOptionsEnd(); 72 73 Units units = honee->units; 74 75 rc = fabs(rc) * units->meter; 76 theta0 *= units->Kelvin; 77 thetaC *= units->Kelvin; 78 P0 *= units->Pascal; 79 N *= (1. / units->second); 80 for (PetscInt i = 0; i < dim; i++) center[i] *= units->meter; 81 82 density_current_ctx->theta0 = theta0; 83 density_current_ctx->thetaC = thetaC; 84 density_current_ctx->P0 = P0; 85 density_current_ctx->N = N; 86 density_current_ctx->rc = rc; 87 PetscCall(PetscArraycpy(density_current_ctx->center, center, 3)); 88 PetscCall(PetscArraycpy(density_current_ctx->dc_axis, dc_axis, 3)); 89 90 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 91 density_current_ctx->newtonian_ctx = *newtonian_ig_ctx; 92 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 93 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &density_current_qfctx)); 94 PetscCallCeed( 95 ceed, CeedQFunctionContextSetData(density_current_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*density_current_ctx), density_current_ctx)); 96 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(density_current_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 97 problem->ics.qfctx = density_current_qfctx; 98 PetscFunctionReturn(PETSC_SUCCESS); 99 } 100