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