1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Miscellaneous utility functions 10 11 #include <ceed.h> 12 #include <petscdm.h> 13 #include <petscsf.h> 14 #include <petscts.h> 15 16 #include "../navierstokes.h" 17 #include "../qfunctions/mass.h" 18 19 PetscErrorCode ICs_FixMultiplicity(DM dm, CeedData ceed_data, User user, Vec Q_loc, Vec Q, CeedScalar time) { 20 Ceed ceed = user->ceed; 21 PetscFunctionBeginUser; 22 23 // --------------------------------------------------------------------------- 24 // Update time for evaluation 25 // --------------------------------------------------------------------------- 26 if (user->phys->ics_time_label) PetscCallCeed(ceed, CeedOperatorSetContextDouble(ceed_data->op_ics_ctx->op, user->phys->ics_time_label, &time)); 27 28 // --------------------------------------------------------------------------- 29 // ICs 30 // --------------------------------------------------------------------------- 31 // -- CEED Restriction 32 CeedVector q0_ceed; 33 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q0_ceed, NULL)); 34 35 // -- Place PETSc vector in CEED vector 36 PetscCall(ApplyCeedOperatorLocalToGlobal(NULL, Q, ceed_data->op_ics_ctx)); 37 38 // --------------------------------------------------------------------------- 39 // Fix multiplicity for output of ICs 40 // --------------------------------------------------------------------------- 41 // -- CEED Restriction 42 CeedVector mult_vec; 43 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &mult_vec, NULL)); 44 45 // -- Place PETSc vector in CEED vector 46 PetscMemType m_mem_type; 47 Vec multiplicity_loc; 48 PetscCall(DMGetLocalVector(dm, &multiplicity_loc)); 49 PetscCall(VecP2C(multiplicity_loc, &m_mem_type, mult_vec)); 50 51 // -- Get multiplicity 52 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, mult_vec)); 53 54 // -- Restore vectors 55 PetscCall(VecC2P(mult_vec, m_mem_type, multiplicity_loc)); 56 57 // -- Local-to-Global 58 Vec multiplicity; 59 PetscCall(DMGetGlobalVector(dm, &multiplicity)); 60 PetscCall(VecZeroEntries(multiplicity)); 61 PetscCall(DMLocalToGlobal(dm, multiplicity_loc, ADD_VALUES, multiplicity)); 62 63 // -- Fix multiplicity 64 PetscCall(VecPointwiseDivide(Q, Q, multiplicity)); 65 PetscCall(VecPointwiseDivide(Q_loc, Q_loc, multiplicity_loc)); 66 67 // -- Restore vectors 68 PetscCall(DMRestoreLocalVector(dm, &multiplicity_loc)); 69 PetscCall(DMRestoreGlobalVector(dm, &multiplicity)); 70 71 // Cleanup 72 PetscCallCeed(ceed, CeedVectorDestroy(&mult_vec)); 73 PetscCallCeed(ceed, CeedVectorDestroy(&q0_ceed)); 74 75 PetscFunctionReturn(PETSC_SUCCESS); 76 } 77 78 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, Vec cell_geom_FVM, 79 Vec grad_FVM) { 80 Vec Qbc, boundary_mask; 81 PetscFunctionBegin; 82 83 // Mask (zero) Strong BC entries 84 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 85 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 86 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 87 88 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 89 PetscCall(VecAXPY(Q_loc, 1., Qbc)); 90 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 91 92 PetscFunctionReturn(PETSC_SUCCESS); 93 } 94 95 // @brief Load vector from binary file, possibly with embedded solution time and step number 96 PetscErrorCode LoadFluidsBinaryVec(MPI_Comm comm, PetscViewer viewer, Vec Q, PetscReal *time, PetscInt *step_number) { 97 PetscInt file_step_number; 98 PetscInt32 token; 99 PetscReal file_time; 100 PetscFunctionBeginUser; 101 102 // Attempt 103 PetscCall(PetscViewerBinaryRead(viewer, &token, 1, NULL, PETSC_INT32)); 104 if (token == FLUIDS_FILE_TOKEN_32 || token == FLUIDS_FILE_TOKEN_64 || 105 token == FLUIDS_FILE_TOKEN) { // New style format; we're reading a file with step number and time in the header 106 PetscCall(PetscViewerBinaryRead(viewer, &file_step_number, 1, NULL, PETSC_INT)); 107 PetscCall(PetscViewerBinaryRead(viewer, &file_time, 1, NULL, PETSC_REAL)); 108 if (time) *time = file_time; 109 if (step_number) *step_number = file_step_number; 110 } else if (token == VEC_FILE_CLASSID) { // Legacy format of just the vector, encoded as [VEC_FILE_CLASSID, length, ] 111 PetscInt length, N; 112 PetscCall(PetscViewerBinaryRead(viewer, &length, 1, NULL, PETSC_INT)); 113 PetscCall(VecGetSize(Q, &N)); 114 PetscCheck(length == N, comm, PETSC_ERR_ARG_INCOMP, "File Vec has length %" PetscInt_FMT " but DM has global Vec size %" PetscInt_FMT, length, N); 115 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE)); 116 } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Not a fluids header token or a PETSc Vec in file"); 117 118 // Load Q from existent solution 119 PetscCall(VecLoad(Q, viewer)); 120 121 PetscFunctionReturn(PETSC_SUCCESS); 122 } 123 124 // Compare reference solution values with current test run for CI 125 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) { 126 Vec Qref; 127 PetscViewer viewer; 128 PetscReal error, Qrefnorm; 129 MPI_Comm comm = PetscObjectComm((PetscObject)Q); 130 PetscFunctionBegin; 131 132 // Read reference file 133 PetscCall(VecDuplicate(Q, &Qref)); 134 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->test_file_path, FILE_MODE_READ, &viewer)); 135 PetscCall(LoadFluidsBinaryVec(comm, viewer, Qref, NULL, NULL)); 136 137 // Compute error with respect to reference solution 138 PetscCall(VecAXPY(Q, -1.0, Qref)); 139 PetscCall(VecNorm(Qref, NORM_MAX, &Qrefnorm)); 140 PetscCall(VecScale(Q, 1. / Qrefnorm)); 141 PetscCall(VecNorm(Q, NORM_MAX, &error)); 142 143 // Check error 144 if (error > app_ctx->test_tol) { 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test failed with error norm %g\n", (double)error)); 146 } 147 148 // Cleanup 149 PetscCall(PetscViewerDestroy(&viewer)); 150 PetscCall(VecDestroy(&Qref)); 151 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154 155 // Get error for problems with exact solutions 156 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, PetscScalar final_time) { 157 PetscInt loc_nodes; 158 Vec Q_exact, Q_exact_loc; 159 PetscReal rel_error, norm_error, norm_exact; 160 PetscFunctionBegin; 161 162 // Get exact solution at final time 163 PetscCall(DMCreateGlobalVector(dm, &Q_exact)); 164 PetscCall(DMGetLocalVector(dm, &Q_exact_loc)); 165 PetscCall(VecGetSize(Q_exact_loc, &loc_nodes)); 166 PetscCall(ICs_FixMultiplicity(dm, ceed_data, user, Q_exact_loc, Q_exact, final_time)); 167 168 // Get |exact solution - obtained solution| 169 PetscCall(VecNorm(Q_exact, NORM_1, &norm_exact)); 170 PetscCall(VecAXPY(Q, -1.0, Q_exact)); 171 PetscCall(VecNorm(Q, NORM_1, &norm_error)); 172 173 // Compute relative error 174 rel_error = norm_error / norm_exact; 175 176 // Output relative error 177 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error: %g\n", (double)rel_error)); 178 // Cleanup 179 PetscCall(DMRestoreLocalVector(dm, &Q_exact_loc)); 180 PetscCall(VecDestroy(&Q_exact)); 181 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 // Post-processing 186 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, ProblemData *problem, User user, Vec Q, PetscScalar final_time) { 187 PetscInt steps; 188 TSConvergedReason reason; 189 PetscFunctionBegin; 190 191 // Print relative error 192 if (problem->non_zero_time && user->app_ctx->test_type == TESTTYPE_NONE) { 193 PetscCall(GetError_NS(ceed_data, dm, user, Q, final_time)); 194 } 195 196 // Print final time and number of steps 197 PetscCall(TSGetStepNumber(ts, &steps)); 198 PetscCall(TSGetConvergedReason(ts, &reason)); 199 if (user->app_ctx->test_type == TESTTYPE_NONE) { 200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time integrator %s on time step %" PetscInt_FMT " with final time %g\n", TSConvergedReasons[reason], 201 steps, (double)final_time)); 202 } 203 204 // Output numerical values from command line 205 PetscCall(VecViewFromOptions(Q, NULL, "-vec_view")); 206 207 // Compare reference solution values with current test run for CI 208 if (user->app_ctx->test_type == TESTTYPE_SOLVER) { 209 PetscCall(RegressionTests_NS(user->app_ctx, Q)); 210 } 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 const PetscInt32 FLUIDS_FILE_TOKEN = 0xceedf00; // for backwards compatibility 215 const PetscInt32 FLUIDS_FILE_TOKEN_32 = 0xceedf32; 216 const PetscInt32 FLUIDS_FILE_TOKEN_64 = 0xceedf64; 217 218 // Gather initial Q values in case of continuation of simulation 219 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q) { 220 PetscViewer viewer; 221 222 PetscFunctionBegin; 223 224 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_file, FILE_MODE_READ, &viewer)); 225 PetscCall(LoadFluidsBinaryVec(comm, viewer, Q, &app_ctx->cont_time, &app_ctx->cont_steps)); 226 PetscCall(PetscViewerDestroy(&viewer)); 227 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 // Record boundary values from initial condition 232 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc) { 233 Vec Qbc, boundary_mask; 234 PetscFunctionBegin; 235 236 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 237 PetscCall(VecCopy(Q_loc, Qbc)); 238 PetscCall(VecZeroEntries(Q_loc)); 239 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, Q_loc)); 240 PetscCall(VecAXPY(Qbc, -1., Q_loc)); 241 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 242 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS)); 243 244 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 245 PetscCall(DMGetGlobalVector(dm, &Q)); 246 PetscCall(VecZeroEntries(boundary_mask)); 247 PetscCall(VecSet(Q, 1.0)); 248 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, boundary_mask)); 249 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 250 251 PetscFunctionReturn(PETSC_SUCCESS); 252 } 253 254 // Free a plain data context that was allocated using PETSc; returning libCEED error codes 255 int FreeContextPetsc(void *data) { 256 if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed"); 257 return CEED_ERROR_SUCCESS; 258 } 259 260 // Return mass qfunction specification for number of components N 261 PetscErrorCode CreateMassQFunction(Ceed ceed, CeedInt N, CeedInt q_data_size, CeedQFunction *qf) { 262 PetscFunctionBeginUser; 263 264 switch (N) { 265 case 1: 266 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Mass_1, Mass_1_loc, qf)); 267 break; 268 case 5: 269 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Mass_5, Mass_5_loc, qf)); 270 break; 271 case 7: 272 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Mass_7, Mass_7_loc, qf)); 273 break; 274 case 9: 275 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Mass_9, Mass_9_loc, qf)); 276 break; 277 case 22: 278 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Mass_22, Mass_22_loc, qf)); 279 break; 280 default: 281 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Could not find mass qfunction of size %d", N); 282 } 283 284 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf, "u", N, CEED_EVAL_INTERP)); 285 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf, "qdata", q_data_size, CEED_EVAL_NONE)); 286 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf, "v", N, CEED_EVAL_INTERP)); 287 PetscFunctionReturn(PETSC_SUCCESS); 288 } 289 290 PetscErrorCode NodalProjectionDataDestroy(NodalProjectionData context) { 291 PetscFunctionBeginUser; 292 if (context == NULL) PetscFunctionReturn(PETSC_SUCCESS); 293 294 PetscCall(DMDestroy(&context->dm)); 295 PetscCall(KSPDestroy(&context->ksp)); 296 297 PetscCall(OperatorApplyContextDestroy(context->l2_rhs_ctx)); 298 299 PetscCall(PetscFree(context)); 300 301 PetscFunctionReturn(PETSC_SUCCESS); 302 } 303 304 /* 305 * @brief Open a PHASTA *.dat file, grabbing dimensions and file pointer 306 * 307 * This function opens the file specified by `path` using `PetscFOpen` and passes the file pointer in `fp`. 308 * It is not closed in this function, thus `fp` must be closed sometime after this function has been called (using `PetscFClose` for example). 309 * 310 * Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space. 311 * 312 * @param[in] comm MPI_Comm for the program 313 * @param[in] path Path to the file 314 * @param[in] char_array_len Length of the character array that should contain each line 315 * @param[out] dims Dimensions of the file, taken from the first line of the file 316 * @param[out] fp File pointer to the opened file 317 */ 318 PetscErrorCode PHASTADatFileOpen(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], const PetscInt char_array_len, PetscInt dims[2], 319 FILE **fp) { 320 int ndims; 321 char line[char_array_len]; 322 char **array; 323 324 PetscFunctionBeginUser; 325 PetscCall(PetscFOpen(comm, path, "r", fp)); 326 PetscCall(PetscSynchronizedFGets(comm, *fp, char_array_len, line)); 327 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 328 PetscCheck(ndims == 2, comm, PETSC_ERR_FILE_UNEXPECTED, "Found %d dimensions instead of 2 on the first line of %s", ndims, path); 329 330 for (PetscInt i = 0; i < ndims; i++) dims[i] = atoi(array[i]); 331 PetscCall(PetscStrToArrayDestroy(ndims, array)); 332 333 PetscFunctionReturn(PETSC_SUCCESS); 334 } 335 336 /* 337 * @brief Get the number of rows for the PHASTA file at path. 338 * 339 * Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space. 340 * 341 * @param[in] comm MPI_Comm for the program 342 * @param[in] path Path to the file 343 * @param[out] nrows Number of rows 344 */ 345 PetscErrorCode PHASTADatFileGetNRows(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscInt *nrows) { 346 const PetscInt char_array_len = 512; 347 PetscInt dims[2]; 348 FILE *fp; 349 350 PetscFunctionBeginUser; 351 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 352 *nrows = dims[0]; 353 PetscCall(PetscFClose(comm, fp)); 354 355 PetscFunctionReturn(PETSC_SUCCESS); 356 } 357 358 PetscErrorCode PHASTADatFileReadToArrayReal(MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscReal array[]) { 359 PetscInt dims[2]; 360 int ndims; 361 FILE *fp; 362 const PetscInt char_array_len = 512; 363 char line[char_array_len]; 364 char **row_array; 365 PetscFunctionBeginUser; 366 367 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 368 369 for (PetscInt i = 0; i < dims[0]; i++) { 370 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 371 PetscCall(PetscStrToArray(line, ' ', &ndims, &row_array)); 372 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED, 373 "Line %" PetscInt_FMT " of %s does not contain enough columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]); 374 375 for (PetscInt j = 0; j < dims[1]; j++) { 376 array[i * dims[1] + j] = (PetscReal)atof(row_array[j]); 377 } 378 } 379 380 PetscCall(PetscFClose(comm, fp)); 381 382 PetscFunctionReturn(PETSC_SUCCESS); 383 } 384 385 PetscLogEvent FLUIDS_CeedOperatorApply; 386 PetscLogEvent FLUIDS_CeedOperatorAssemble; 387 PetscLogEvent FLUIDS_CeedOperatorAssembleDiagonal; 388 PetscLogEvent FLUIDS_CeedOperatorAssemblePointBlockDiagonal; 389 static PetscClassId libCEED_classid; 390 391 PetscErrorCode RegisterLogEvents() { 392 PetscFunctionBeginUser; 393 PetscCall(PetscClassIdRegister("libCEED", &libCEED_classid)); 394 PetscCall(PetscLogEventRegister("CeedOpApply", libCEED_classid, &FLUIDS_CeedOperatorApply)); 395 PetscCall(PetscLogEventRegister("CeedOpAsm", libCEED_classid, &FLUIDS_CeedOperatorAssemble)); 396 PetscCall(PetscLogEventRegister("CeedOpAsmD", libCEED_classid, &FLUIDS_CeedOperatorAssembleDiagonal)); 397 PetscCall(PetscLogEventRegister("CeedOpAsmPBD", libCEED_classid, &FLUIDS_CeedOperatorAssemblePointBlockDiagonal)); 398 PetscFunctionReturn(PETSC_SUCCESS); 399 } 400 401 /** 402 @brief Translate array of CeedInt to PetscInt. 403 If the types differ, `array_ceed` is freed with `free()` and `array_petsc` is allocated with `malloc()`. 404 Caller is responsible for freeing `array_petsc` with `free()`. 405 406 @param[in] num_entries Number of array entries 407 @param[in,out] array_ceed Array of CeedInts 408 @param[out] array_petsc Array of PetscInts 409 **/ 410 PetscErrorCode IntArrayC2P(PetscInt num_entries, CeedInt **array_ceed, PetscInt **array_petsc) { 411 CeedInt int_c = 0; 412 PetscInt int_p = 0; 413 414 PetscFunctionBeginUser; 415 if (sizeof(int_c) == sizeof(int_p)) { 416 *array_petsc = (PetscInt *)*array_ceed; 417 } else { 418 *array_petsc = malloc(num_entries * sizeof(PetscInt)); 419 for (PetscInt i = 0; i < num_entries; i++) (*array_petsc)[i] = (*array_ceed)[i]; 420 free(*array_ceed); 421 } 422 *array_ceed = NULL; 423 424 PetscFunctionReturn(PETSC_SUCCESS); 425 } 426 427 /** 428 @brief Translate array of PetscInt to CeedInt. 429 If the types differ, `array_petsc` is freed with `PetscFree()` and `array_ceed` is allocated with `PetscMalloc1()`. 430 Caller is responsible for freeing `array_ceed` with `PetscFree()`. 431 432 @param[in] num_entries Number of array entries 433 @param[in,out] array_petsc Array of PetscInts 434 @param[out] array_ceed Array of CeedInts 435 **/ 436 PetscErrorCode IntArrayP2C(PetscInt num_entries, PetscInt **array_petsc, CeedInt **array_ceed) { 437 CeedInt int_c = 0; 438 PetscInt int_p = 0; 439 440 PetscFunctionBeginUser; 441 if (sizeof(int_c) == sizeof(int_p)) { 442 *array_ceed = (CeedInt *)*array_petsc; 443 } else { 444 PetscCall(PetscMalloc1(num_entries, array_ceed)); 445 for (PetscInt i = 0; i < num_entries; i++) (*array_ceed)[i] = (*array_petsc)[i]; 446 PetscCall(PetscFree(*array_petsc)); 447 } 448 *array_petsc = NULL; 449 450 PetscFunctionReturn(PETSC_SUCCESS); 451 } 452 453 // Print information about the given simulation run 454 PetscErrorCode PrintRunInfo(User user, Physics phys_ctx, ProblemData *problem, MPI_Comm comm) { 455 Ceed ceed = user->ceed; 456 PetscFunctionBeginUser; 457 // Header and rank 458 char host_name[PETSC_MAX_PATH_LEN]; 459 PetscMPIInt rank, comm_size; 460 PetscCall(PetscGetHostName(host_name, sizeof host_name)); 461 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 462 PetscCallMPI(MPI_Comm_size(comm, &comm_size)); 463 PetscCall(PetscPrintf(comm, 464 "\n-- Navier-Stokes solver - libCEED + PETSc --\n" 465 " MPI:\n" 466 " Host Name : %s\n" 467 " Total ranks : %d\n", 468 host_name, comm_size)); 469 470 // Problem specific info 471 PetscCall(problem->print_info(user, problem, user->app_ctx)); 472 473 // libCEED 474 const char *used_resource; 475 CeedMemType mem_type_backend; 476 PetscCallCeed(ceed, CeedGetResource(user->ceed, &used_resource)); 477 PetscCallCeed(ceed, CeedGetPreferredMemType(user->ceed, &mem_type_backend)); 478 PetscCall(PetscPrintf(comm, 479 " libCEED:\n" 480 " libCEED Backend : %s\n" 481 " libCEED Backend MemType : %s\n", 482 used_resource, CeedMemTypes[mem_type_backend])); 483 // PETSc 484 char box_faces_str[PETSC_MAX_PATH_LEN] = "3,3,3"; 485 if (problem->dim == 2) box_faces_str[3] = '\0'; 486 PetscCall(PetscOptionsGetString(NULL, NULL, "-dm_plex_box_faces", box_faces_str, sizeof(box_faces_str), NULL)); 487 MatType mat_type; 488 VecType vec_type; 489 PetscCall(DMGetMatType(user->dm, &mat_type)); 490 PetscCall(DMGetVecType(user->dm, &vec_type)); 491 PetscCall(PetscPrintf(comm, 492 " PETSc:\n" 493 " Box Faces : %s\n" 494 " DM MatType : %s\n" 495 " DM VecType : %s\n" 496 " Time Stepping Scheme : %s\n", 497 box_faces_str, mat_type, vec_type, phys_ctx->implicit ? "implicit" : "explicit")); 498 if (user->app_ctx->cont_steps) { 499 PetscCall(PetscPrintf(comm, 500 " Continue:\n" 501 " Filename: : %s\n" 502 " Step: : %" PetscInt_FMT "\n" 503 " Time: : %g\n", 504 user->app_ctx->cont_file, user->app_ctx->cont_steps, user->app_ctx->cont_time)); 505 } 506 // Mesh 507 const PetscInt num_comp_q = 5; 508 PetscInt glob_dofs, owned_dofs, local_dofs; 509 const CeedInt num_P = user->app_ctx->degree + 1, num_Q = num_P + user->app_ctx->q_extra; 510 // -- Get global size 511 PetscCall(DMGetGlobalVectorInfo(user->dm, &owned_dofs, &glob_dofs, NULL)); 512 // -- Get local size 513 PetscCall(DMGetLocalVectorInfo(user->dm, &local_dofs, NULL, NULL)); 514 PetscCall(PetscPrintf(comm, 515 " Mesh:\n" 516 " Number of 1D Basis Nodes (P) : %" CeedInt_FMT "\n" 517 " Number of 1D Quadrature Points (Q) : %" CeedInt_FMT "\n" 518 " Global DoFs : %" PetscInt_FMT "\n" 519 " DoFs per node : %" PetscInt_FMT "\n" 520 " Global %" PetscInt_FMT "-DoF nodes : %" PetscInt_FMT "\n", 521 num_P, num_Q, glob_dofs, num_comp_q, num_comp_q, glob_dofs / num_comp_q)); 522 // -- Get Partition Statistics 523 PetscCall(PetscPrintf(comm, " Partition: (min,max,median,max/median)\n")); 524 { 525 PetscInt *gather_buffer = NULL; 526 PetscInt part_owned_dofs[3], part_local_dofs[3], part_boundary_dofs[3], part_neighbors[3]; 527 PetscInt median_index = comm_size % 2 ? comm_size / 2 : comm_size / 2 - 1; 528 if (!rank) PetscCall(PetscMalloc1(comm_size, &gather_buffer)); 529 530 PetscCallMPI(MPI_Gather(&owned_dofs, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 531 if (!rank) { 532 PetscCall(PetscSortInt(comm_size, gather_buffer)); 533 part_owned_dofs[0] = gather_buffer[0]; // min 534 part_owned_dofs[1] = gather_buffer[comm_size - 1]; // max 535 part_owned_dofs[2] = gather_buffer[median_index]; // median 536 PetscReal part_owned_dof_ratio = (PetscReal)part_owned_dofs[1] / (PetscReal)part_owned_dofs[2]; 537 PetscCall(PetscPrintf( 538 comm, " Global Vector %" PetscInt_FMT "-DoF nodes : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", num_comp_q, 539 part_owned_dofs[0] / num_comp_q, part_owned_dofs[1] / num_comp_q, part_owned_dofs[2] / num_comp_q, part_owned_dof_ratio)); 540 } 541 542 PetscCallMPI(MPI_Gather(&local_dofs, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 543 if (!rank) { 544 PetscCall(PetscSortInt(comm_size, gather_buffer)); 545 part_local_dofs[0] = gather_buffer[0]; // min 546 part_local_dofs[1] = gather_buffer[comm_size - 1]; // max 547 part_local_dofs[2] = gather_buffer[median_index]; // median 548 PetscReal part_local_dof_ratio = (PetscReal)part_local_dofs[1] / (PetscReal)part_local_dofs[2]; 549 PetscCall(PetscPrintf( 550 comm, " Local Vector %" PetscInt_FMT "-DoF nodes : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", num_comp_q, 551 part_local_dofs[0] / num_comp_q, part_local_dofs[1] / num_comp_q, part_local_dofs[2] / num_comp_q, part_local_dof_ratio)); 552 } 553 554 if (comm_size != 1) { 555 PetscInt num_remote_roots_total = 0, num_remote_leaves_total = 0, num_ghost_interface_ranks = 0, num_owned_interface_ranks = 0; 556 { 557 PetscSF sf; 558 PetscInt nrranks, niranks; 559 const PetscInt *roffset, *rmine, *rremote, *ioffset, *irootloc; 560 const PetscMPIInt *rranks, *iranks; 561 PetscCall(DMGetSectionSF(user->dm, &sf)); 562 PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote)); 563 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc)); 564 for (PetscInt i = 0; i < nrranks; i++) { 565 if (rranks[i] == rank) continue; // Ignore same-part global->local transfers 566 num_remote_roots_total += roffset[i + 1] - roffset[i]; 567 num_ghost_interface_ranks++; 568 } 569 for (PetscInt i = 0; i < niranks; i++) { 570 if (iranks[i] == rank) continue; 571 num_remote_leaves_total += ioffset[i + 1] - ioffset[i]; 572 num_owned_interface_ranks++; 573 } 574 } 575 PetscCallMPI(MPI_Gather(&num_remote_roots_total, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 576 if (!rank) { 577 PetscCall(PetscSortInt(comm_size, gather_buffer)); 578 part_boundary_dofs[0] = gather_buffer[0]; // min 579 part_boundary_dofs[1] = gather_buffer[comm_size - 1]; // max 580 part_boundary_dofs[2] = gather_buffer[median_index]; // median 581 PetscReal part_shared_dof_ratio = (PetscReal)part_boundary_dofs[1] / (PetscReal)part_boundary_dofs[2]; 582 PetscCall(PetscPrintf( 583 comm, " Ghost Interface %" PetscInt_FMT "-DoF nodes : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", 584 num_comp_q, part_boundary_dofs[0] / num_comp_q, part_boundary_dofs[1] / num_comp_q, part_boundary_dofs[2] / num_comp_q, 585 part_shared_dof_ratio)); 586 } 587 588 PetscCallMPI(MPI_Gather(&num_ghost_interface_ranks, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 589 if (!rank) { 590 PetscCall(PetscSortInt(comm_size, gather_buffer)); 591 part_neighbors[0] = gather_buffer[0]; // min 592 part_neighbors[1] = gather_buffer[comm_size - 1]; // max 593 part_neighbors[2] = gather_buffer[median_index]; // median 594 PetscReal part_neighbors_ratio = (PetscReal)part_neighbors[1] / (PetscReal)part_neighbors[2]; 595 PetscCall(PetscPrintf(comm, " Ghost Interface Ranks : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", 596 part_neighbors[0], part_neighbors[1], part_neighbors[2], part_neighbors_ratio)); 597 } 598 599 PetscCallMPI(MPI_Gather(&num_remote_leaves_total, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 600 if (!rank) { 601 PetscCall(PetscSortInt(comm_size, gather_buffer)); 602 part_boundary_dofs[0] = gather_buffer[0]; // min 603 part_boundary_dofs[1] = gather_buffer[comm_size - 1]; // max 604 part_boundary_dofs[2] = gather_buffer[median_index]; // median 605 PetscReal part_shared_dof_ratio = (PetscReal)part_boundary_dofs[1] / (PetscReal)part_boundary_dofs[2]; 606 PetscCall(PetscPrintf( 607 comm, " Owned Interface %" PetscInt_FMT "-DoF nodes : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", 608 num_comp_q, part_boundary_dofs[0] / num_comp_q, part_boundary_dofs[1] / num_comp_q, part_boundary_dofs[2] / num_comp_q, 609 part_shared_dof_ratio)); 610 } 611 612 PetscCallMPI(MPI_Gather(&num_owned_interface_ranks, 1, MPIU_INT, gather_buffer, 1, MPIU_INT, 0, comm)); 613 if (!rank) { 614 PetscCall(PetscSortInt(comm_size, gather_buffer)); 615 part_neighbors[0] = gather_buffer[0]; // min 616 part_neighbors[1] = gather_buffer[comm_size - 1]; // max 617 part_neighbors[2] = gather_buffer[median_index]; // median 618 PetscReal part_neighbors_ratio = (PetscReal)part_neighbors[1] / (PetscReal)part_neighbors[2]; 619 PetscCall(PetscPrintf(comm, " Owned Interface Ranks : %" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ", %f\n", 620 part_neighbors[0], part_neighbors[1], part_neighbors[2], part_neighbors_ratio)); 621 } 622 } 623 624 if (!rank) PetscCall(PetscFree(gather_buffer)); 625 } 626 627 PetscFunctionReturn(PETSC_SUCCESS); 628 } 629