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 SHOCKTUBE 6 7 #include "../qfunctions/shocktube.h" 8 9 #include <ceed.h> 10 #include <petscdm.h> 11 12 #include <navierstokes.h> 13 14 static PetscErrorCode PRINT_SHOCKTUBE(Honee honee, ProblemData problem, AppCtx app_ctx) { 15 MPI_Comm comm = honee->comm; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscPrintf(comm, 19 " Problem:\n" 20 " Problem Name : %s\n", 21 app_ctx->problem_name)); 22 PetscFunctionReturn(PETSC_SUCCESS); 23 } 24 25 static const char *const component_names[] = {"Density", "MomentumX", "MomentumY", "MomentumZ", "TotalEnergy"}; 26 27 PetscErrorCode NS_SHOCKTUBE(ProblemData problem, DM dm, void *ctx) { 28 SetupContextShock setup_context; 29 Honee honee = *(Honee *)ctx; 30 MPI_Comm comm = honee->comm; 31 Ceed ceed = honee->ceed; 32 PetscBool implicit; 33 PetscBool yzb; 34 PetscInt stab; 35 ShockTubeContext shocktube_ctx; 36 CeedQFunctionContext shocktube_qfctx, ics_qfctx; 37 38 PetscFunctionBeginUser; 39 PetscCall(PetscNew(&setup_context)); 40 PetscCall(PetscNew(&shocktube_ctx)); 41 42 // ------------------------------------------------------ 43 // Create the QFunction context 44 // ------------------------------------------------------ 45 // Driver section initial conditions 46 CeedScalar P_high = 1.0; // Pa 47 CeedScalar rho_high = 1.0; // kg/m^3 48 // Driven section initial conditions 49 CeedScalar P_low = 0.1; // Pa 50 CeedScalar rho_low = 0.125; // kg/m^3 51 // Stabilization parameter 52 CeedScalar c_tau = 0.5; // -, based on Hughes et al (2010) 53 // Tuning parameters for the YZB shock capturing 54 CeedScalar Cyzb = 0.1; // -, used in approximation of (Na),x 55 CeedScalar Byzb = 2.0; // -, 1 for smooth shocks 56 // 2 for sharp shocks 57 PetscReal domain_min[3], domain_max[3], domain_size[3]; 58 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 59 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 60 61 // ------------------------------------------------------ 62 // Command line Options 63 // ------------------------------------------------------ 64 PetscOptionsBegin(comm, NULL, "Options for SHOCKTUBE problem", NULL); 65 66 // -- Numerical formulation options 67 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL)); 68 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL)); 69 PetscCall(PetscOptionsScalar("-c_tau", "Stabilization constant", NULL, c_tau, &c_tau, NULL)); 70 PetscCall(PetscOptionsBool("-yzb", "Use YZB discontinuity capturing", NULL, yzb = PETSC_FALSE, &yzb, NULL)); 71 72 // -- Warnings 73 if (stab == STAB_SUPG) { 74 PetscCall(PetscPrintf(comm, "Warning! -stab supg not implemented for the shocktube problem. \n")); 75 } 76 if (yzb && implicit) { 77 PetscCall(PetscPrintf(comm, "Warning! -yzb only implemented for explicit timestepping. \n")); 78 } 79 80 PetscOptionsEnd(); 81 82 // ------------------------------------------------------ 83 // Set up the QFunction context 84 // ------------------------------------------------------ 85 // -- Scale variables to desired units 86 CeedScalar mid_point = 0.5 * (domain_size[0] + domain_min[0]); 87 88 // -- Setup Context 89 setup_context->mid_point = mid_point; 90 setup_context->time = 0.0; 91 setup_context->P_high = P_high; 92 setup_context->rho_high = rho_high; 93 setup_context->P_low = P_low; 94 setup_context->rho_low = rho_low; 95 96 // -- QFunction Context 97 honee->phys->implicit = implicit; 98 shocktube_ctx->implicit = implicit; 99 shocktube_ctx->stabilization = stab; 100 shocktube_ctx->yzb = yzb; 101 shocktube_ctx->Cyzb = Cyzb; 102 shocktube_ctx->Byzb = Byzb; 103 shocktube_ctx->c_tau = c_tau; 104 105 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &ics_qfctx)); 106 PetscCallCeed(ceed, CeedQFunctionContextSetData(ics_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*setup_context), setup_context)); 107 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(ics_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 108 109 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &shocktube_qfctx)); 110 PetscCallCeed(ceed, CeedQFunctionContextSetData(shocktube_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*shocktube_ctx), shocktube_ctx)); 111 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(shocktube_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 112 113 // ------------------------------------------------------ 114 // SET UP SHOCKTUBE 115 // ------------------------------------------------------ 116 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsShockTube, .qf_loc = ICsShockTube_loc, .qfctx = ics_qfctx}; 117 problem->apply_vol_rhs = (HoneeQFSpec){.qf_func_ptr = EulerShockTube, .qf_loc = EulerShockTube_loc, .qfctx = shocktube_qfctx}; 118 problem->num_comps_jac_data = 0; 119 problem->compute_exact_solution_error = PETSC_FALSE; 120 problem->print_info = PRINT_SHOCKTUBE; 121 122 problem->num_components = 5; 123 PetscCall(PetscMalloc1(problem->num_components, &problem->component_names)); 124 for (PetscInt i = 0; i < 5; i++) PetscCall(PetscStrallocpy(component_names[i], &problem->component_names[i])); 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127