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