1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 /// @file
9 /// Utility functions for setting up DENSITY_CURRENT
10
11 #include "../qfunctions/densitycurrent.h"
12
13 #include <ceed.h>
14 #include <petscdm.h>
15
16 #include "../navierstokes.h"
17
NS_DENSITY_CURRENT(ProblemData problem,DM dm,void * ctx,SimpleBC bc)18 PetscErrorCode NS_DENSITY_CURRENT(ProblemData problem, DM dm, void *ctx, SimpleBC bc) {
19 User user = *(User *)ctx;
20 MPI_Comm comm = user->comm;
21 Ceed ceed = user->ceed;
22 DensityCurrentContext dc_ctx;
23 CeedQFunctionContext density_current_context;
24 NewtonianIdealGasContext newtonian_ig_ctx;
25
26 PetscFunctionBeginUser;
27 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx, bc));
28 PetscCall(PetscCalloc1(1, &dc_ctx));
29 // ------------------------------------------------------
30 // SET UP DENSITY_CURRENT
31 // ------------------------------------------------------
32 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
33 problem->ics.qfunction = ICsDC;
34 problem->ics.qfunction_loc = ICsDC_loc;
35
36 // ------------------------------------------------------
37 // Create the libCEED context
38 // ------------------------------------------------------
39 CeedScalar theta0 = 300.; // K
40 CeedScalar thetaC = -15.; // K
41 CeedScalar P0 = 1.e5; // Pa
42 CeedScalar N = 0.01; // 1/s
43 CeedScalar rc = 1000.; // m (Radius of bubble)
44 PetscReal center[3], dc_axis[3] = {0, 0, 0};
45 PetscReal domain_min[3], domain_max[3], domain_size[3];
46 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
47 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
48
49 // ------------------------------------------------------
50 // Command line Options
51 // ------------------------------------------------------
52 PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem", NULL);
53 PetscCall(PetscOptionsScalar("-theta0", "Reference potential temperature", NULL, theta0, &theta0, NULL));
54 PetscCall(PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", NULL, thetaC, &thetaC, NULL));
55 PetscCall(PetscOptionsScalar("-P0", "Atmospheric pressure", NULL, P0, &P0, NULL));
56 PetscCall(PetscOptionsScalar("-N", "Brunt-Vaisala frequency", NULL, N, &N, NULL));
57 PetscCall(PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", NULL, rc, &rc, NULL));
58 for (PetscInt i = 0; i < 3; i++) center[i] = .5 * domain_size[i];
59 PetscInt n = problem->dim;
60 PetscCall(PetscOptionsRealArray("-center", "Location of bubble center", NULL, center, &n, NULL));
61 n = problem->dim;
62 PetscCall(PetscOptionsRealArray("-dc_axis",
63 "Axis of density current cylindrical anomaly, "
64 "or {0,0,0} for spherically symmetric",
65 NULL, dc_axis, &n, NULL));
66 {
67 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + PetscSqr(dc_axis[1]) + PetscSqr(dc_axis[2]));
68 if (norm > 0) {
69 for (PetscInt i = 0; i < 3; i++) dc_axis[i] /= norm;
70 }
71 }
72
73 PetscOptionsEnd();
74
75 PetscScalar meter = user->units->meter;
76 PetscScalar second = user->units->second;
77 PetscScalar Kelvin = user->units->Kelvin;
78 PetscScalar Pascal = user->units->Pascal;
79 rc = fabs(rc) * meter;
80 theta0 *= Kelvin;
81 thetaC *= Kelvin;
82 P0 *= Pascal;
83 N *= (1. / second);
84 for (PetscInt i = 0; i < 3; i++) center[i] *= meter;
85
86 dc_ctx->theta0 = theta0;
87 dc_ctx->thetaC = thetaC;
88 dc_ctx->P0 = P0;
89 dc_ctx->N = N;
90 dc_ctx->rc = rc;
91 dc_ctx->center[0] = center[0];
92 dc_ctx->center[1] = center[1];
93 dc_ctx->center[2] = center[2];
94 dc_ctx->dc_axis[0] = dc_axis[0];
95 dc_ctx->dc_axis[1] = dc_axis[1];
96 dc_ctx->dc_axis[2] = dc_axis[2];
97
98 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
99 dc_ctx->newtonian_ctx = *newtonian_ig_ctx;
100 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx));
101 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &density_current_context));
102 PetscCallCeed(ceed, CeedQFunctionContextSetData(density_current_context, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*dc_ctx), dc_ctx));
103 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(density_current_context, CEED_MEM_HOST, FreeContextPetsc));
104 problem->ics.qfunction_context = density_current_context;
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107