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