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