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