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 "../navierstokes.h" 12 13 PetscErrorCode ICs_FixMultiplicity(DM dm, CeedData ceed_data, User user, Vec Q_loc, Vec Q, CeedScalar time) { 14 PetscFunctionBeginUser; 15 16 // --------------------------------------------------------------------------- 17 // Update time for evaluation 18 // --------------------------------------------------------------------------- 19 if (user->phys->ics_time_label) CeedOperatorContextSetDouble(ceed_data->op_ics, user->phys->ics_time_label, &time); 20 21 // --------------------------------------------------------------------------- 22 // ICs 23 // --------------------------------------------------------------------------- 24 // -- CEED Restriction 25 CeedVector q0_ceed; 26 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q0_ceed, NULL); 27 28 // -- Place PETSc vector in CEED vector 29 CeedScalar *q0; 30 PetscMemType q0_mem_type; 31 PetscCall(VecGetArrayAndMemType(Q_loc, (PetscScalar **)&q0, &q0_mem_type)); 32 CeedVectorSetArray(q0_ceed, MemTypeP2C(q0_mem_type), CEED_USE_POINTER, q0); 33 34 // -- Apply CEED Operator 35 CeedOperatorApply(ceed_data->op_ics, ceed_data->x_coord, q0_ceed, CEED_REQUEST_IMMEDIATE); 36 37 // -- Restore vectors 38 CeedVectorTakeArray(q0_ceed, MemTypeP2C(q0_mem_type), NULL); 39 PetscCall(VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q0)); 40 41 // -- Local-to-Global 42 PetscCall(VecZeroEntries(Q)); 43 PetscCall(DMLocalToGlobal(dm, Q_loc, ADD_VALUES, Q)); 44 45 // --------------------------------------------------------------------------- 46 // Fix multiplicity for output of ICs 47 // --------------------------------------------------------------------------- 48 // -- CEED Restriction 49 CeedVector mult_vec; 50 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &mult_vec, NULL); 51 52 // -- Place PETSc vector in CEED vector 53 CeedScalar *mult; 54 PetscMemType m_mem_type; 55 Vec multiplicity_loc; 56 PetscCall(DMGetLocalVector(dm, &multiplicity_loc)); 57 PetscCall(VecGetArrayAndMemType(multiplicity_loc, (PetscScalar **)&mult, &m_mem_type)); 58 CeedVectorSetArray(mult_vec, MemTypeP2C(m_mem_type), CEED_USE_POINTER, mult); 59 60 // -- Get multiplicity 61 CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, mult_vec); 62 63 // -- Restore vectors 64 CeedVectorTakeArray(mult_vec, MemTypeP2C(m_mem_type), NULL); 65 PetscCall(VecRestoreArrayReadAndMemType(multiplicity_loc, (const PetscScalar **)&mult)); 66 67 // -- Local-to-Global 68 Vec multiplicity; 69 PetscCall(DMGetGlobalVector(dm, &multiplicity)); 70 PetscCall(VecZeroEntries(multiplicity)); 71 PetscCall(DMLocalToGlobal(dm, multiplicity_loc, ADD_VALUES, multiplicity)); 72 73 // -- Fix multiplicity 74 PetscCall(VecPointwiseDivide(Q, Q, multiplicity)); 75 PetscCall(VecPointwiseDivide(Q_loc, Q_loc, multiplicity_loc)); 76 77 // -- Restore vectors 78 PetscCall(DMRestoreLocalVector(dm, &multiplicity_loc)); 79 PetscCall(DMRestoreGlobalVector(dm, &multiplicity)); 80 81 // Cleanup 82 CeedVectorDestroy(&mult_vec); 83 CeedVectorDestroy(&q0_ceed); 84 85 PetscFunctionReturn(0); 86 } 87 88 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, Vec cell_geom_FVM, 89 Vec grad_FVM) { 90 Vec Qbc, boundary_mask; 91 PetscFunctionBegin; 92 93 // Mask (zero) Dirichlet entries 94 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 95 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 96 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 97 98 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 99 PetscCall(VecAXPY(Q_loc, 1., Qbc)); 100 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 101 102 PetscFunctionReturn(0); 103 } 104 105 // Compare reference solution values with current test run for CI 106 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) { 107 Vec Qref; 108 PetscViewer viewer; 109 PetscReal error, Qrefnorm; 110 PetscFunctionBegin; 111 112 // Read reference file 113 PetscCall(VecDuplicate(Q, &Qref)); 114 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)Q), app_ctx->file_path, FILE_MODE_READ, &viewer)); 115 PetscCall(VecLoad(Qref, viewer)); 116 117 // Compute error with respect to reference solution 118 PetscCall(VecAXPY(Q, -1.0, Qref)); 119 PetscCall(VecNorm(Qref, NORM_MAX, &Qrefnorm)); 120 PetscCall(VecScale(Q, 1. / Qrefnorm)); 121 PetscCall(VecNorm(Q, NORM_MAX, &error)); 122 123 // Check error 124 if (error > app_ctx->test_tol) { 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test failed with error norm %g\n", (double)error)); 126 } 127 128 // Cleanup 129 PetscCall(PetscViewerDestroy(&viewer)); 130 PetscCall(VecDestroy(&Qref)); 131 132 PetscFunctionReturn(0); 133 } 134 135 // Get error for problems with exact solutions 136 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, PetscScalar final_time) { 137 PetscInt loc_nodes; 138 Vec Q_exact, Q_exact_loc; 139 PetscReal rel_error, norm_error, norm_exact; 140 PetscFunctionBegin; 141 142 // Get exact solution at final time 143 PetscCall(DMCreateGlobalVector(dm, &Q_exact)); 144 PetscCall(DMGetLocalVector(dm, &Q_exact_loc)); 145 PetscCall(VecGetSize(Q_exact_loc, &loc_nodes)); 146 PetscCall(ICs_FixMultiplicity(dm, ceed_data, user, Q_exact_loc, Q_exact, final_time)); 147 148 // Get |exact solution - obtained solution| 149 PetscCall(VecNorm(Q_exact, NORM_1, &norm_exact)); 150 PetscCall(VecAXPY(Q, -1.0, Q_exact)); 151 PetscCall(VecNorm(Q, NORM_1, &norm_error)); 152 153 // Compute relative error 154 rel_error = norm_error / norm_exact; 155 156 // Output relative error 157 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error: %g\n", (double)rel_error)); 158 // Cleanup 159 PetscCall(DMRestoreLocalVector(dm, &Q_exact_loc)); 160 PetscCall(VecDestroy(&Q_exact)); 161 162 PetscFunctionReturn(0); 163 } 164 165 // Post-processing 166 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, ProblemData *problem, User user, Vec Q, PetscScalar final_time) { 167 PetscInt steps; 168 TSConvergedReason reason; 169 PetscFunctionBegin; 170 171 // Print relative error 172 if (problem->non_zero_time && !user->app_ctx->test_mode) { 173 PetscCall(GetError_NS(ceed_data, dm, user, Q, final_time)); 174 } 175 176 // Print final time and number of steps 177 PetscCall(TSGetStepNumber(ts, &steps)); 178 PetscCall(TSGetConvergedReason(ts, &reason)); 179 if (!user->app_ctx->test_mode) { 180 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time integrator %s on time step %" PetscInt_FMT " with final time %g\n", TSConvergedReasons[reason], 181 steps, (double)final_time)); 182 } 183 184 // Output numerical values from command line 185 PetscCall(VecViewFromOptions(Q, NULL, "-vec_view")); 186 187 // Compare reference solution values with current test run for CI 188 if (user->app_ctx->test_mode) { 189 PetscCall(RegressionTests_NS(user->app_ctx, Q)); 190 } 191 192 PetscFunctionReturn(0); 193 } 194 195 const PetscInt FLUIDS_FILE_TOKEN = 0xceedf00; 196 197 // Gather initial Q values in case of continuation of simulation 198 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q) { 199 PetscViewer viewer; 200 PetscInt token, step_number; 201 PetscReal time; 202 203 PetscFunctionBegin; 204 205 // Read input 206 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_file, FILE_MODE_READ, &viewer)); 207 208 // Attempt 209 PetscCall(PetscViewerBinaryRead(viewer, &token, 1, NULL, PETSC_INT)); 210 if (token == FLUIDS_FILE_TOKEN) { // New style format; we're reading a file with step number and time in the header 211 PetscCall(PetscViewerBinaryRead(viewer, &step_number, 1, NULL, PETSC_INT)); 212 PetscCall(PetscViewerBinaryRead(viewer, &time, 1, NULL, PETSC_REAL)); 213 app_ctx->cont_steps = step_number; 214 app_ctx->cont_time = time; 215 } else if (token == VEC_FILE_CLASSID) { // Legacy format of just the vector, encoded as [VEC_FILE_CLASSID, length, ] 216 PetscInt length, N; 217 PetscCall(PetscViewerBinaryRead(viewer, &length, 1, NULL, PETSC_INT)); 218 PetscCall(VecGetSize(Q, &N)); 219 PetscCheck(length == N, comm, PETSC_ERR_ARG_INCOMP, "File Vec has length %" PetscInt_FMT " but DM has global Vec size %" PetscInt_FMT, length, N); 220 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE)); 221 } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Not a fluids header token or a PETSc Vec in file"); 222 223 // Load Q from existent solution 224 PetscCall(VecLoad(Q, viewer)); 225 226 // Cleanup 227 PetscCall(PetscViewerDestroy(&viewer)); 228 229 PetscFunctionReturn(0); 230 } 231 232 // Record boundary values from initial condition 233 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc) { 234 Vec Qbc, boundary_mask; 235 PetscFunctionBegin; 236 237 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 238 PetscCall(VecCopy(Q_loc, Qbc)); 239 PetscCall(VecZeroEntries(Q_loc)); 240 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, Q_loc)); 241 PetscCall(VecAXPY(Qbc, -1., Q_loc)); 242 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 243 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS)); 244 245 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 246 PetscCall(DMGetGlobalVector(dm, &Q)); 247 PetscCall(VecZeroEntries(boundary_mask)); 248 PetscCall(VecSet(Q, 1.0)); 249 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, boundary_mask)); 250 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 251 252 PetscFunctionReturn(0); 253 } 254 255 // Free a plain data context that was allocated using PETSc; returning libCEED error codes 256 int FreeContextPetsc(void *data) { 257 if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed"); 258 return CEED_ERROR_SUCCESS; 259 } 260