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 SetupContext setup_context; 19 User user = *(User *)ctx; 20 MPI_Comm comm = PETSC_COMM_WORLD; 21 22 PetscFunctionBeginUser; 23 ierr = NS_NEWTONIAN_IG(problem, dm, ctx); CHKERRQ(ierr); 24 // ------------------------------------------------------ 25 // SET UP DENSITY_CURRENT 26 // ------------------------------------------------------ 27 problem->ics.qfunction = ICsDC; 28 problem->ics.qfunction_loc = ICsDC_loc; 29 problem->bc = Exact_DC; 30 CeedQFunctionContextGetData(problem->ics.qfunction_context, CEED_MEM_HOST, 31 &setup_context); 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 ierr = DMGetBoundingBox(dm, domain_min, domain_max); 44 CHKERRQ(ierr); 45 for (PetscInt i = 0; i < 3; i++) 46 domain_size[i] = domain_max[i] - domain_min[i]; 47 48 // ------------------------------------------------------ 49 // Command line Options 50 // ------------------------------------------------------ 51 PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem", NULL); 52 ierr = PetscOptionsScalar("-theta0", "Reference potential temperature", NULL, 53 theta0, &theta0, NULL); 54 CHKERRQ(ierr); 55 ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature", 56 NULL, thetaC, &thetaC, NULL); 57 CHKERRQ(ierr); 58 ierr = PetscOptionsScalar("-P0", "Atmospheric pressure", NULL, P0, &P0, NULL); 59 CHKERRQ(ierr); 60 ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency", NULL, N, &N, NULL); 61 CHKERRQ(ierr); 62 ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble", 63 NULL, rc, &rc, NULL); 64 CHKERRQ(ierr); 65 for (PetscInt i = 0; i < 3; i++) 66 center[i] = .5 * domain_size[i]; 67 PetscInt n = problem->dim; 68 ierr = PetscOptionsRealArray("-center", "Location of bubble center", NULL, 69 center, &n, NULL); 70 CHKERRQ(ierr); 71 n = problem->dim; 72 ierr = PetscOptionsRealArray("-dc_axis", 73 "Axis of density current cylindrical anomaly, " 74 "or {0,0,0} for spherically symmetric", 75 NULL, dc_axis, &n, NULL); 76 CHKERRQ(ierr); 77 { 78 PetscReal norm = PetscSqrtReal(PetscSqr(dc_axis[0]) + PetscSqr(dc_axis[1]) + 79 PetscSqr(dc_axis[2])); 80 if (norm > 0) { 81 for (PetscInt i = 0; i < 3; i++) 82 dc_axis[i] /= norm; 83 } 84 } 85 86 PetscOptionsEnd(); 87 88 PetscScalar meter = user->units->meter; 89 PetscScalar second = user->units->second; 90 PetscScalar Kelvin = user->units->Kelvin; 91 PetscScalar Pascal = user->units->Pascal; 92 rc = fabs(rc) * meter; 93 theta0 *= Kelvin; 94 thetaC *= Kelvin; 95 P0 *= Pascal; 96 N *= (1. / second); 97 for (PetscInt i = 0; i < 3; i++) 98 center[i] *= meter; 99 100 setup_context->theta0 = theta0; 101 setup_context->thetaC = thetaC; 102 setup_context->P0 = P0; 103 setup_context->N = N; 104 setup_context->rc = rc; 105 setup_context->center[0] = center[0]; 106 setup_context->center[1] = center[1]; 107 setup_context->center[2] = center[2]; 108 setup_context->dc_axis[0] = dc_axis[0]; 109 setup_context->dc_axis[1] = dc_axis[1]; 110 setup_context->dc_axis[2] = dc_axis[2]; 111 112 problem->bc_ctx = 113 setup_context; // This is bad, context data should only be accessed via Get/Restore 114 CeedQFunctionContextRestoreData(problem->ics.qfunction_context, &setup_context); 115 116 PetscFunctionReturn(0); 117 } 118