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