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, SimpleBC bc) { 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, bc)); 25 PetscCall(PetscCalloc1(1, &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 libCEED 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 PetscScalar meter = honee->units->meter; 74 PetscScalar second = honee->units->second; 75 PetscScalar Kelvin = honee->units->Kelvin; 76 PetscScalar Pascal = honee->units->Pascal; 77 rc = fabs(rc) * meter; 78 theta0 *= Kelvin; 79 thetaC *= Kelvin; 80 P0 *= Pascal; 81 N *= (1. / second); 82 for (PetscInt i = 0; i < dim; i++) center[i] *= meter; 83 84 density_current_ctx->theta0 = theta0; 85 density_current_ctx->thetaC = thetaC; 86 density_current_ctx->P0 = P0; 87 density_current_ctx->N = N; 88 density_current_ctx->rc = rc; 89 PetscCall(PetscArraycpy(density_current_ctx->center, center, 3)); 90 PetscCall(PetscArraycpy(density_current_ctx->dc_axis, dc_axis, 3)); 91 92 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 93 density_current_ctx->newtonian_ctx = *newtonian_ig_ctx; 94 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 95 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &density_current_qfctx)); 96 PetscCallCeed( 97 ceed, CeedQFunctionContextSetData(density_current_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*density_current_ctx), density_current_ctx)); 98 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(density_current_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 99 problem->ics.qfctx = density_current_qfctx; 100 PetscFunctionReturn(PETSC_SUCCESS); 101 } 102