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
NS_DENSITY_CURRENT(ProblemData problem,DM dm,void * ctx)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->newt_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(ceed, CeedQFunctionContextSetData(density_current_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*density_current_ctx),
89 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