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; 37 38 PetscFunctionBeginUser; 39 PetscCall(PetscCalloc1(1, &setup_context)); 40 PetscCall(PetscCalloc1(1, &shocktube_ctx)); 41 42 // ------------------------------------------------------ 43 // SET UP SHOCKTUBE 44 // ------------------------------------------------------ 45 problem->ics.qf_func_ptr = ICsShockTube; 46 problem->ics.qf_loc = ICsShockTube_loc; 47 problem->apply_vol_rhs.qf_func_ptr = EulerShockTube; 48 problem->apply_vol_rhs.qf_loc = EulerShockTube_loc; 49 problem->apply_vol_ifunction.qf_func_ptr = NULL; 50 problem->apply_vol_ifunction.qf_loc = NULL; 51 problem->num_comps_jac_data = 0; 52 problem->compute_exact_solution_error = PETSC_FALSE; 53 problem->print_info = PRINT_SHOCKTUBE; 54 55 problem->num_components = 5; 56 PetscCall(PetscMalloc1(problem->num_components, &problem->component_names)); 57 for (PetscInt i = 0; i < 5; i++) PetscCall(PetscStrallocpy(component_names[i], &problem->component_names[i])); 58 59 // ------------------------------------------------------ 60 // Create the QFunction context 61 // ------------------------------------------------------ 62 // Driver section initial conditions 63 CeedScalar P_high = 1.0; // Pa 64 CeedScalar rho_high = 1.0; // kg/m^3 65 // Driven section initial conditions 66 CeedScalar P_low = 0.1; // Pa 67 CeedScalar rho_low = 0.125; // kg/m^3 68 // Stabilization parameter 69 CeedScalar c_tau = 0.5; // -, based on Hughes et al (2010) 70 // Tuning parameters for the YZB shock capturing 71 CeedScalar Cyzb = 0.1; // -, used in approximation of (Na),x 72 CeedScalar Byzb = 2.0; // -, 1 for smooth shocks 73 // 2 for sharp shocks 74 PetscReal domain_min[3], domain_max[3], domain_size[3]; 75 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 76 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 77 78 // ------------------------------------------------------ 79 // Command line Options 80 // ------------------------------------------------------ 81 PetscOptionsBegin(comm, NULL, "Options for SHOCKTUBE problem", NULL); 82 83 // -- Numerical formulation options 84 PetscCall(PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", NULL, implicit = PETSC_FALSE, &implicit, NULL)); 85 PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL)); 86 PetscCall(PetscOptionsScalar("-c_tau", "Stabilization constant", NULL, c_tau, &c_tau, NULL)); 87 PetscCall(PetscOptionsBool("-yzb", "Use YZB discontinuity capturing", NULL, yzb = PETSC_FALSE, &yzb, NULL)); 88 89 // -- Warnings 90 if (stab == STAB_SUPG) { 91 PetscCall(PetscPrintf(comm, "Warning! -stab supg not implemented for the shocktube problem. \n")); 92 } 93 if (yzb && implicit) { 94 PetscCall(PetscPrintf(comm, "Warning! -yzb only implemented for explicit timestepping. \n")); 95 } 96 97 PetscOptionsEnd(); 98 99 // ------------------------------------------------------ 100 // Set up the QFunction context 101 // ------------------------------------------------------ 102 // -- Scale variables to desired units 103 Units units = honee->units; 104 for (PetscInt i = 0; i < 3; i++) { 105 domain_size[i] *= units->meter; 106 domain_min[i] *= units->meter; 107 } 108 CeedScalar mid_point = 0.5 * (domain_size[0] + domain_min[0]); 109 110 // -- Setup Context 111 setup_context->mid_point = mid_point; 112 setup_context->time = 0.0; 113 setup_context->P_high = P_high; 114 setup_context->rho_high = rho_high; 115 setup_context->P_low = P_low; 116 setup_context->rho_low = rho_low; 117 118 // -- QFunction Context 119 honee->phys->implicit = implicit; 120 shocktube_ctx->implicit = implicit; 121 shocktube_ctx->stabilization = stab; 122 shocktube_ctx->yzb = yzb; 123 shocktube_ctx->Cyzb = Cyzb; 124 shocktube_ctx->Byzb = Byzb; 125 shocktube_ctx->c_tau = c_tau; 126 127 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &problem->ics.qfctx)); 128 PetscCallCeed(ceed, CeedQFunctionContextSetData(problem->ics.qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*setup_context), setup_context)); 129 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(problem->ics.qfctx, CEED_MEM_HOST, FreeContextPetsc)); 130 131 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &shocktube_qfctx)); 132 PetscCallCeed(ceed, CeedQFunctionContextSetData(shocktube_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*shocktube_ctx), shocktube_ctx)); 133 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(shocktube_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 134 problem->apply_vol_rhs.qfctx = shocktube_qfctx; 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137