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