1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #ifndef libceed_fluids_examples_navier_stokes_h 18 #define libceed_fluids_examples_navier_stokes_h 19 20 #include <ceed.h> 21 #include <petscdm.h> 22 #include <petscdmplex.h> 23 #include <petscsys.h> 24 #include <petscts.h> 25 #include <stdbool.h> 26 27 // ----------------------------------------------------------------------------- 28 // PETSc Version 29 // ----------------------------------------------------------------------------- 30 #if PETSC_VERSION_LT(3,17,0) 31 #error "PETSc v3.17 or later is required" 32 #endif 33 34 // ----------------------------------------------------------------------------- 35 // Enums 36 // ----------------------------------------------------------------------------- 37 // Translate PetscMemType to CeedMemType 38 static inline CeedMemType MemTypeP2C(PetscMemType mem_type) { 39 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 40 } 41 42 // Advection - Wind Options 43 typedef enum { 44 WIND_ROTATION = 0, 45 WIND_TRANSLATION = 1, 46 } WindType; 47 static const char *const WindTypes[] = { 48 "rotation", 49 "translation", 50 "WindType", "WIND_", NULL 51 }; 52 53 // Advection - Bubble Types 54 typedef enum { 55 BUBBLE_SPHERE = 0, // dim=3 56 BUBBLE_CYLINDER = 1, // dim=2 57 } BubbleType; 58 static const char *const BubbleTypes[] = { 59 "sphere", 60 "cylinder", 61 "BubbleType", "BUBBLE_", NULL 62 }; 63 64 // Advection - Bubble Continuity Types 65 typedef enum { 66 BUBBLE_CONTINUITY_SMOOTH = 0, // Original continuous, smooth shape 67 BUBBLE_CONTINUITY_BACK_SHARP = 1, // Discontinuous, sharp back half shape 68 BUBBLE_CONTINUITY_THICK = 2, // Define a finite thickness 69 } BubbleContinuityType; 70 static const char *const BubbleContinuityTypes[] = { 71 "smooth", 72 "back_sharp", 73 "thick", 74 "BubbleContinuityType", "BUBBLE_CONTINUITY_", NULL 75 }; 76 77 // Euler - test cases 78 typedef enum { 79 EULER_TEST_ISENTROPIC_VORTEX = 0, 80 EULER_TEST_1 = 1, 81 EULER_TEST_2 = 2, 82 EULER_TEST_3 = 3, 83 EULER_TEST_4 = 4, 84 EULER_TEST_5 = 5, 85 } EulerTestType; 86 static const char *const EulerTestTypes[] = { 87 "isentropic_vortex", 88 "test_1", 89 "test_2", 90 "test_3", 91 "test_4", 92 "test_5", 93 "EulerTestType", "EULER_TEST_", NULL 94 }; 95 96 // Stabilization methods 97 typedef enum { 98 STAB_NONE = 0, 99 STAB_SU = 1, // Streamline Upwind 100 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 101 } StabilizationType; 102 static const char *const StabilizationTypes[] = { 103 "none", 104 "SU", 105 "SUPG", 106 "StabilizationType", "STAB_", NULL 107 }; 108 109 // ----------------------------------------------------------------------------- 110 // Structs 111 // ----------------------------------------------------------------------------- 112 // Structs declarations 113 typedef struct AppCtx_private *AppCtx; 114 typedef struct CeedData_private *CeedData; 115 typedef struct User_private *User; 116 typedef struct Units_private *Units; 117 typedef struct SimpleBC_private *SimpleBC; 118 typedef struct Physics_private *Physics; 119 120 // Application context from user command line options 121 struct AppCtx_private { 122 // libCEED arguments 123 char ceed_resource[PETSC_MAX_PATH_LEN]; // libCEED backend 124 PetscInt degree; 125 PetscInt q_extra; 126 // Post-processing arguments 127 PetscInt output_freq; 128 PetscInt viz_refine; 129 PetscInt cont_steps; 130 char output_dir[PETSC_MAX_PATH_LEN]; 131 // Problem type arguments 132 PetscFunctionList problems; 133 char problem_name[PETSC_MAX_PATH_LEN]; 134 // Test mode arguments 135 PetscBool test_mode; 136 PetscScalar test_tol; 137 char file_path[PETSC_MAX_PATH_LEN]; 138 }; 139 140 // libCEED data struct 141 struct CeedData_private { 142 CeedVector x_coord, q_data; 143 CeedQFunctionContext setup_context, newt_ig_context, advection_context, 144 euler_context; 145 CeedQFunction qf_setup_vol, qf_ics, qf_rhs_vol, qf_ifunction_vol, 146 qf_setup_sur, qf_apply_inflow, qf_apply_outflow; 147 CeedBasis basis_x, basis_xc, basis_q, basis_x_sur, basis_q_sur; 148 CeedElemRestriction elem_restr_x, elem_restr_q, elem_restr_qd_i; 149 CeedOperator op_setup_vol, op_ics; 150 }; 151 152 // PETSc user data 153 struct User_private { 154 MPI_Comm comm; 155 DM dm; 156 DM dm_viz; 157 Mat interp_viz; 158 Ceed ceed; 159 Units units; 160 Vec M; 161 Physics phys; 162 AppCtx app_ctx; 163 CeedVector q_ceed, q_dot_ceed, g_ceed; 164 CeedOperator op_rhs_vol, op_rhs, op_ifunction_vol, op_ifunction; 165 }; 166 167 // Units 168 struct Units_private { 169 // fundamental units 170 PetscScalar meter; 171 PetscScalar kilogram; 172 PetscScalar second; 173 PetscScalar Kelvin; 174 // derived units 175 PetscScalar Pascal; 176 PetscScalar J_per_kg_K; 177 PetscScalar m_per_squared_s; 178 PetscScalar W_per_m_K; 179 PetscScalar Joule; 180 }; 181 182 // Boundary conditions 183 struct SimpleBC_private { 184 PetscInt num_wall, // Number of faces with wall BCs 185 wall_comps[5], // An array of constrained component numbers 186 num_comps, 187 num_slip[3], // Number of faces with slip BCs 188 num_inflow, 189 num_outflow; 190 PetscInt walls[16], slips[3][16], inflows[16], outflows[16]; 191 PetscBool user_bc; 192 }; 193 194 // Initial conditions 195 #ifndef setup_context_struct 196 #define setup_context_struct 197 typedef struct SetupContext_ *SetupContext; 198 struct SetupContext_ { 199 CeedScalar theta0; 200 CeedScalar thetaC; 201 CeedScalar P0; 202 CeedScalar N; 203 CeedScalar cv; 204 CeedScalar cp; 205 CeedScalar g; 206 CeedScalar rc; 207 CeedScalar lx; 208 CeedScalar ly; 209 CeedScalar lz; 210 CeedScalar center[3]; 211 CeedScalar dc_axis[3]; 212 CeedScalar wind[3]; 213 CeedScalar time; 214 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 215 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 216 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 217 }; 218 #endif 219 220 // DENSITY_CURRENT 221 #ifndef dc_context_struct 222 #define dc_context_struct 223 typedef struct DCContext_ *DCContext; 224 struct DCContext_ { 225 CeedScalar lambda; 226 CeedScalar mu; 227 CeedScalar k; 228 CeedScalar cv; 229 CeedScalar cp; 230 CeedScalar g; 231 CeedScalar c_tau; 232 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 233 }; 234 #endif 235 236 // EULER_VORTEX 237 #ifndef euler_context_struct 238 #define euler_context_struct 239 typedef struct EulerContext_ *EulerContext; 240 struct EulerContext_ { 241 CeedScalar center[3]; 242 CeedScalar curr_time; 243 CeedScalar vortex_strength; 244 CeedScalar c_tau; 245 CeedScalar mean_velocity[3]; 246 bool implicit; 247 int euler_test; 248 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 249 }; 250 #endif 251 252 // ADVECTION and ADVECTION2D 253 #ifndef advection_context_struct 254 #define advection_context_struct 255 typedef struct AdvectionContext_ *AdvectionContext; 256 struct AdvectionContext_ { 257 CeedScalar CtauS; 258 CeedScalar strong_form; 259 CeedScalar E_wind; 260 bool implicit; 261 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 262 }; 263 #endif 264 265 // Newtonian Ideal Gas 266 #ifndef newtonian_context_struct 267 #define newtonian_context_struct 268 typedef struct NewtonianIdealGasContext_ *NewtonianIdealGasContext; 269 struct NewtonianIdealGasContext_ { 270 CeedScalar lambda; 271 CeedScalar mu; 272 CeedScalar k; 273 CeedScalar cv; 274 CeedScalar cp; 275 CeedScalar g; 276 CeedScalar c_tau; 277 StabilizationType stabilization; 278 }; 279 #endif 280 281 // Struct that contains all enums and structs used for the physics of all problems 282 struct Physics_private { 283 NewtonianIdealGasContext newtonian_ig_ctx; 284 EulerContext euler_ctx; 285 AdvectionContext advection_ctx; 286 WindType wind_type; 287 BubbleType bubble_type; 288 BubbleContinuityType bubble_continuity_type; 289 EulerTestType euler_test; 290 StabilizationType stab; 291 PetscBool implicit; 292 PetscBool has_curr_time; 293 PetscBool has_neumann; 294 }; 295 296 // Problem specific data 297 // *INDENT-OFF* 298 typedef struct { 299 CeedInt dim, q_data_size_vol, q_data_size_sur; 300 CeedScalar dm_scale; 301 CeedQFunctionUser setup_vol, setup_sur, ics, apply_vol_rhs, apply_vol_ifunction, 302 apply_inflow, apply_outflow; 303 const char *setup_vol_loc, *setup_sur_loc, *ics_loc, 304 *apply_vol_rhs_loc, *apply_vol_ifunction_loc, *apply_inflow_loc, *apply_outflow_loc; 305 bool non_zero_time; 306 PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, 307 PetscScalar[], void *); 308 PetscErrorCode (*setup_ctx)(Ceed, CeedData, AppCtx, SetupContext, Physics); 309 PetscErrorCode (*print_info)(Physics, SetupContext, AppCtx); 310 } ProblemData; 311 // *INDENT-ON* 312 313 // ----------------------------------------------------------------------------- 314 // Set up problems 315 // ----------------------------------------------------------------------------- 316 // Set up function for each problem 317 extern PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, 318 void *setup_ctx, void *ctx); 319 extern PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, 320 void *setup_ctx, void *ctx); 321 extern PetscErrorCode NS_EULER_VORTEX(ProblemData *problem, DM dm, 322 void *setup_ctx, void *ctx); 323 extern PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *setup_ctx, 324 void *ctx); 325 extern PetscErrorCode NS_ADVECTION2D(ProblemData *problem, DM dm, 326 void *setup_ctx, void *ctx); 327 328 // Set up context for each problem 329 extern PetscErrorCode SetupContext_NEWTONIAN_IG(Ceed ceed, CeedData ceed_data, 330 AppCtx app_ctx, SetupContext setup_ctx, Physics phys); 331 332 extern PetscErrorCode SetupContext_DENSITY_CURRENT(Ceed ceed, 333 CeedData ceed_data, AppCtx app_ctx, SetupContext setup_ctx, Physics phys); 334 335 extern PetscErrorCode SetupContext_EULER_VORTEX(Ceed ceed, CeedData ceed_data, 336 AppCtx app_ctx, SetupContext setup_ctx, Physics phys); 337 338 extern PetscErrorCode SetupContext_ADVECTION(Ceed ceed, CeedData ceed_data, 339 AppCtx app_ctx, SetupContext setup_ctx, Physics phys); 340 341 extern PetscErrorCode SetupContext_ADVECTION2D(Ceed ceed, CeedData ceed_data, 342 AppCtx app_ctx, SetupContext setup_ctx, Physics phys); 343 344 // Boundary condition function for each problem 345 extern PetscErrorCode BC_DENSITY_CURRENT(DM dm, SimpleBC bc, Physics phys, 346 void *setup_ctx); 347 348 extern PetscErrorCode BC_EULER_VORTEX(DM dm, SimpleBC bc, Physics phys, 349 void *setup_ctx); 350 351 extern PetscErrorCode BC_ADVECTION(DM dm, SimpleBC bc, Physics phys, 352 void *setup_ctx); 353 354 extern PetscErrorCode BC_ADVECTION2D(DM dm, SimpleBC bc, Physics phys, 355 void *setup_ctx); 356 357 // Print function for each problem 358 extern PetscErrorCode PRINT_DENSITY_CURRENT(Physics phys, 359 SetupContext setup_ctx, AppCtx app_ctx); 360 361 extern PetscErrorCode PRINT_EULER_VORTEX(Physics phys, SetupContext setup_ctx, 362 AppCtx app_ctx); 363 364 extern PetscErrorCode PRINT_ADVECTION(Physics phys, SetupContext setup_ctx, 365 AppCtx app_ctx); 366 367 extern PetscErrorCode PRINT_ADVECTION2D(Physics phys, SetupContext setup_ctx, 368 AppCtx app_ctx); 369 370 // ----------------------------------------------------------------------------- 371 // libCEED functions 372 // ----------------------------------------------------------------------------- 373 // Utility function - essential BC dofs are encoded in closure indices as -(i+1). 374 PetscInt Involute(PetscInt i); 375 376 // Utility function to create local CEED restriction 377 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, 378 DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr); 379 380 // Utility function to get Ceed Restriction for each domain 381 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 382 DMLabel domain_label, PetscInt value, 383 CeedInt Q, CeedInt q_data_size, 384 CeedElemRestriction *elem_restr_q, 385 CeedElemRestriction *elem_restr_x, 386 CeedElemRestriction *elem_restr_qd_i); 387 388 // Utility function to create CEED Composite Operator for the entire domain 389 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, 390 CeedData ceed_data, Physics phys, 391 CeedOperator op_apply_vol, CeedInt height, 392 CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur, 393 CeedOperator *op_apply); 394 395 PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, 396 AppCtx app_ctx, ProblemData *problem, SimpleBC bc); 397 398 // ----------------------------------------------------------------------------- 399 // Time-stepping functions 400 // ----------------------------------------------------------------------------- 401 // Compute mass matrix for explicit scheme 402 PetscErrorCode ComputeLumpedMassMatrix(Ceed ceed, DM dm, CeedData ceed_data, 403 Vec M); 404 405 // RHS (Explicit time-stepper) function setup 406 PetscErrorCode RHS_NS(TS ts, PetscReal t, Vec Q, Vec G, void *user_data); 407 408 // Implicit time-stepper function setup 409 PetscErrorCode IFunction_NS(TS ts, PetscReal t, Vec Q, Vec Q_dot, Vec G, 410 void *user_data); 411 412 // User provided TS Monitor 413 PetscErrorCode TSMonitor_NS(TS ts, PetscInt step_no, PetscReal time, Vec Q, 414 void *ctx); 415 416 // TS: Create, setup, and solve 417 PetscErrorCode TSSolve_NS(DM dm, User user, AppCtx app_ctx, Physics phys, 418 Vec *Q, PetscScalar *f_time, TS *ts); 419 420 // ----------------------------------------------------------------------------- 421 // Setup DM 422 // ----------------------------------------------------------------------------- 423 // Create mesh 424 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, DM *dm); 425 426 // Set up DM 427 PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, 428 SimpleBC bc, Physics phys, void *setup_ctx); 429 430 // Refine DM for high-order viz 431 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, 432 SimpleBC bc, Physics phys, void *setup_ctx); 433 434 // ----------------------------------------------------------------------------- 435 // Process command line options 436 // ----------------------------------------------------------------------------- 437 // Register problems to be available on the command line 438 PetscErrorCode RegisterProblems_NS(AppCtx app_ctx); 439 440 // Process general command line options 441 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx, 442 SimpleBC bc); 443 444 // ----------------------------------------------------------------------------- 445 // Miscellaneous utility functions 446 // ----------------------------------------------------------------------------- 447 PetscErrorCode ICs_FixMultiplicity(DM dm, CeedData ceed_data, Vec Q_loc, Vec Q, 448 CeedScalar time); 449 450 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 451 PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 452 Vec cell_geom_FVM, Vec grad_FVM); 453 454 // Compare reference solution values with current test run for CI 455 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q); 456 457 // Get error for problems with exact solutions 458 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, AppCtx app_ctx, Vec Q, 459 PetscScalar final_time); 460 461 // Post-processing 462 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, 463 ProblemData *problem, AppCtx app_ctx, 464 Vec Q, PetscScalar final_time); 465 466 // -- Gather initial Q values in case of continuation of simulation 467 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q); 468 469 // Record boundary values from initial condition 470 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc); 471 472 // ----------------------------------------------------------------------------- 473 474 #endif // libceed_fluids_examples_navier_stokes_h 475