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 #include "../qfunctions/mass.h" 13 14 PetscErrorCode ICs_FixMultiplicity(DM dm, CeedData ceed_data, User user, Vec Q_loc, Vec Q, CeedScalar time) { 15 PetscFunctionBeginUser; 16 17 // --------------------------------------------------------------------------- 18 // Update time for evaluation 19 // --------------------------------------------------------------------------- 20 if (user->phys->ics_time_label) CeedOperatorContextSetDouble(ceed_data->op_ics, user->phys->ics_time_label, &time); 21 22 // --------------------------------------------------------------------------- 23 // ICs 24 // --------------------------------------------------------------------------- 25 // -- CEED Restriction 26 CeedVector q0_ceed; 27 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q0_ceed, NULL); 28 29 // -- Place PETSc vector in CEED vector 30 PetscMemType q0_mem_type; 31 PetscCall(VecP2C(Q_loc, &q0_mem_type, q0_ceed)); 32 33 // -- Apply CEED Operator 34 CeedOperatorApply(ceed_data->op_ics, ceed_data->x_coord, q0_ceed, CEED_REQUEST_IMMEDIATE); 35 36 // -- Restore vectors 37 PetscCall(VecC2P(q0_ceed, q0_mem_type, Q_loc)); 38 39 // -- Local-to-Global 40 PetscCall(VecZeroEntries(Q)); 41 PetscCall(DMLocalToGlobal(dm, Q_loc, ADD_VALUES, Q)); 42 43 // --------------------------------------------------------------------------- 44 // Fix multiplicity for output of ICs 45 // --------------------------------------------------------------------------- 46 // -- CEED Restriction 47 CeedVector mult_vec; 48 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &mult_vec, NULL); 49 50 // -- Place PETSc vector in CEED vector 51 PetscMemType m_mem_type; 52 Vec multiplicity_loc; 53 PetscCall(DMGetLocalVector(dm, &multiplicity_loc)); 54 PetscCall(VecP2C(multiplicity_loc, &m_mem_type, mult_vec)); 55 56 // -- Get multiplicity 57 CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, mult_vec); 58 59 // -- Restore vectors 60 PetscCall(VecC2P(mult_vec, m_mem_type, multiplicity_loc)); 61 62 // -- Local-to-Global 63 Vec multiplicity; 64 PetscCall(DMGetGlobalVector(dm, &multiplicity)); 65 PetscCall(VecZeroEntries(multiplicity)); 66 PetscCall(DMLocalToGlobal(dm, multiplicity_loc, ADD_VALUES, multiplicity)); 67 68 // -- Fix multiplicity 69 PetscCall(VecPointwiseDivide(Q, Q, multiplicity)); 70 PetscCall(VecPointwiseDivide(Q_loc, Q_loc, multiplicity_loc)); 71 72 // -- Restore vectors 73 PetscCall(DMRestoreLocalVector(dm, &multiplicity_loc)); 74 PetscCall(DMRestoreGlobalVector(dm, &multiplicity)); 75 76 // Cleanup 77 CeedVectorDestroy(&mult_vec); 78 CeedVectorDestroy(&q0_ceed); 79 80 PetscFunctionReturn(0); 81 } 82 83 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, Vec cell_geom_FVM, 84 Vec grad_FVM) { 85 Vec Qbc, boundary_mask; 86 PetscFunctionBegin; 87 88 // Mask (zero) Dirichlet entries 89 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 90 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 91 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 92 93 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 94 PetscCall(VecAXPY(Q_loc, 1., Qbc)); 95 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 96 97 PetscFunctionReturn(0); 98 } 99 100 // Compare reference solution values with current test run for CI 101 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) { 102 Vec Qref; 103 PetscViewer viewer; 104 PetscReal error, Qrefnorm; 105 PetscFunctionBegin; 106 107 // Read reference file 108 PetscCall(VecDuplicate(Q, &Qref)); 109 PetscCall(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)Q), app_ctx->test_file_path, FILE_MODE_READ, &viewer)); 110 PetscCall(VecLoad(Qref, viewer)); 111 112 // Compute error with respect to reference solution 113 PetscCall(VecAXPY(Q, -1.0, Qref)); 114 PetscCall(VecNorm(Qref, NORM_MAX, &Qrefnorm)); 115 PetscCall(VecScale(Q, 1. / Qrefnorm)); 116 PetscCall(VecNorm(Q, NORM_MAX, &error)); 117 118 // Check error 119 if (error > app_ctx->test_tol) { 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test failed with error norm %g\n", (double)error)); 121 } 122 123 // Cleanup 124 PetscCall(PetscViewerDestroy(&viewer)); 125 PetscCall(VecDestroy(&Qref)); 126 127 PetscFunctionReturn(0); 128 } 129 130 // Get error for problems with exact solutions 131 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, PetscScalar final_time) { 132 PetscInt loc_nodes; 133 Vec Q_exact, Q_exact_loc; 134 PetscReal rel_error, norm_error, norm_exact; 135 PetscFunctionBegin; 136 137 // Get exact solution at final time 138 PetscCall(DMCreateGlobalVector(dm, &Q_exact)); 139 PetscCall(DMGetLocalVector(dm, &Q_exact_loc)); 140 PetscCall(VecGetSize(Q_exact_loc, &loc_nodes)); 141 PetscCall(ICs_FixMultiplicity(dm, ceed_data, user, Q_exact_loc, Q_exact, final_time)); 142 143 // Get |exact solution - obtained solution| 144 PetscCall(VecNorm(Q_exact, NORM_1, &norm_exact)); 145 PetscCall(VecAXPY(Q, -1.0, Q_exact)); 146 PetscCall(VecNorm(Q, NORM_1, &norm_error)); 147 148 // Compute relative error 149 rel_error = norm_error / norm_exact; 150 151 // Output relative error 152 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error: %g\n", (double)rel_error)); 153 // Cleanup 154 PetscCall(DMRestoreLocalVector(dm, &Q_exact_loc)); 155 PetscCall(VecDestroy(&Q_exact)); 156 157 PetscFunctionReturn(0); 158 } 159 160 // Post-processing 161 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, ProblemData *problem, User user, Vec Q, PetscScalar final_time) { 162 PetscInt steps; 163 TSConvergedReason reason; 164 PetscFunctionBegin; 165 166 // Print relative error 167 if (problem->non_zero_time && user->app_ctx->test_type == TESTTYPE_NONE) { 168 PetscCall(GetError_NS(ceed_data, dm, user, Q, final_time)); 169 } 170 171 // Print final time and number of steps 172 PetscCall(TSGetStepNumber(ts, &steps)); 173 PetscCall(TSGetConvergedReason(ts, &reason)); 174 if (user->app_ctx->test_type == TESTTYPE_NONE) { 175 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time integrator %s on time step %" PetscInt_FMT " with final time %g\n", TSConvergedReasons[reason], 176 steps, (double)final_time)); 177 } 178 179 // Output numerical values from command line 180 PetscCall(VecViewFromOptions(Q, NULL, "-vec_view")); 181 182 // Compare reference solution values with current test run for CI 183 if (user->app_ctx->test_type == TESTTYPE_SOLVER) { 184 PetscCall(RegressionTests_NS(user->app_ctx, Q)); 185 } 186 187 PetscFunctionReturn(0); 188 } 189 190 const PetscInt FLUIDS_FILE_TOKEN = 0xceedf00; 191 192 // Gather initial Q values in case of continuation of simulation 193 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q) { 194 PetscViewer viewer; 195 PetscInt token, step_number; 196 PetscReal time; 197 198 PetscFunctionBegin; 199 200 // Read input 201 PetscCall(PetscViewerBinaryOpen(comm, app_ctx->cont_file, FILE_MODE_READ, &viewer)); 202 203 // Attempt 204 PetscCall(PetscViewerBinaryRead(viewer, &token, 1, NULL, PETSC_INT)); 205 if (token == FLUIDS_FILE_TOKEN) { // New style format; we're reading a file with step number and time in the header 206 PetscCall(PetscViewerBinaryRead(viewer, &step_number, 1, NULL, PETSC_INT)); 207 PetscCall(PetscViewerBinaryRead(viewer, &time, 1, NULL, PETSC_REAL)); 208 app_ctx->cont_steps = step_number; 209 app_ctx->cont_time = time; 210 } else if (token == VEC_FILE_CLASSID) { // Legacy format of just the vector, encoded as [VEC_FILE_CLASSID, length, ] 211 PetscInt length, N; 212 PetscCall(PetscViewerBinaryRead(viewer, &length, 1, NULL, PETSC_INT)); 213 PetscCall(VecGetSize(Q, &N)); 214 PetscCheck(length == N, comm, PETSC_ERR_ARG_INCOMP, "File Vec has length %" PetscInt_FMT " but DM has global Vec size %" PetscInt_FMT, length, N); 215 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE)); 216 } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Not a fluids header token or a PETSc Vec in file"); 217 218 // Load Q from existent solution 219 PetscCall(VecLoad(Q, viewer)); 220 221 // Cleanup 222 PetscCall(PetscViewerDestroy(&viewer)); 223 224 PetscFunctionReturn(0); 225 } 226 227 // Record boundary values from initial condition 228 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc) { 229 Vec Qbc, boundary_mask; 230 PetscFunctionBegin; 231 232 PetscCall(DMGetNamedLocalVector(dm, "Qbc", &Qbc)); 233 PetscCall(VecCopy(Q_loc, Qbc)); 234 PetscCall(VecZeroEntries(Q_loc)); 235 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, Q_loc)); 236 PetscCall(VecAXPY(Qbc, -1., Q_loc)); 237 PetscCall(DMRestoreNamedLocalVector(dm, "Qbc", &Qbc)); 238 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS)); 239 240 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 241 PetscCall(DMGetGlobalVector(dm, &Q)); 242 PetscCall(VecZeroEntries(boundary_mask)); 243 PetscCall(VecSet(Q, 1.0)); 244 PetscCall(DMGlobalToLocal(dm, Q, INSERT_VALUES, boundary_mask)); 245 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 246 247 PetscFunctionReturn(0); 248 } 249 250 // Free a plain data context that was allocated using PETSc; returning libCEED error codes 251 int FreeContextPetsc(void *data) { 252 if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed"); 253 return CEED_ERROR_SUCCESS; 254 } 255 256 // Return mass qfunction specification for number of components N 257 PetscErrorCode CreateMassQFunction(Ceed ceed, CeedInt N, CeedInt q_data_size, CeedQFunction *qf) { 258 CeedQFunctionUser qfunction_ptr; 259 const char *qfunction_loc; 260 PetscFunctionBeginUser; 261 262 switch (N) { 263 case 1: 264 qfunction_ptr = Mass_1; 265 qfunction_loc = Mass_1_loc; 266 break; 267 case 5: 268 qfunction_ptr = Mass_5; 269 qfunction_loc = Mass_5_loc; 270 break; 271 case 9: 272 qfunction_ptr = Mass_9; 273 qfunction_loc = Mass_9_loc; 274 break; 275 case 22: 276 qfunction_ptr = Mass_22; 277 qfunction_loc = Mass_22_loc; 278 break; 279 default: 280 SETERRQ(PETSC_COMM_WORLD, -1, "Could not find mass qfunction of size %d", N); 281 } 282 CeedQFunctionCreateInterior(ceed, 1, qfunction_ptr, qfunction_loc, qf); 283 284 CeedQFunctionAddInput(*qf, "u", N, CEED_EVAL_INTERP); 285 CeedQFunctionAddInput(*qf, "qdata", q_data_size, CEED_EVAL_NONE); 286 CeedQFunctionAddOutput(*qf, "v", N, CEED_EVAL_INTERP); 287 PetscFunctionReturn(0); 288 } 289 290 /* @brief L^2 Projection of a source FEM function to a target FEM space 291 * 292 * To solve system using a lumped mass matrix, pass a KSP object with ksp_type=preonly, pc_type=jacobi, pc_jacobi_type=rowsum. 293 * 294 * @param[in] source_vec Global Vec of the source FEM function. NULL indicates using rhs_matop_ctx->X_loc 295 * @param[out] target_vec Global Vec of the target (result) FEM function. NULL indicates using rhs_matop_ctx->Y_loc 296 * @param[in] rhs_matop_ctx MatopApplyContext for performing the RHS evaluation 297 * @param[in] ksp KSP for solving the consistent projection problem 298 */ 299 PetscErrorCode ComputeL2Projection(Vec source_vec, Vec target_vec, MatopApplyContext rhs_matop_ctx, KSP ksp) { 300 PetscFunctionBeginUser; 301 302 PetscCall(ApplyLocal_Ceed(source_vec, target_vec, rhs_matop_ctx)); 303 PetscCall(KSPSolve(ksp, target_vec, target_vec)); 304 305 PetscFunctionReturn(0); 306 } 307