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 102 // Note: The BCs must be inserted *before* the other values are inserted into Q_loc 103 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 104 PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 105 Vec cell_geom_FVM, Vec grad_FVM) { 106 107 Vec Qbc; 108 PetscErrorCode ierr; 109 PetscFunctionBegin; 110 111 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 112 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 113 ierr = VecAXPY(Q_loc, 1., Qbc); CHKERRQ(ierr); 114 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 115 116 PetscFunctionReturn(0); 117 } 118 119 // Compare reference solution values with current test run for CI 120 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) { 121 122 Vec Qref; 123 PetscViewer viewer; 124 PetscReal error, Qrefnorm; 125 PetscErrorCode ierr; 126 PetscFunctionBegin; 127 128 // Read reference file 129 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 130 ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)Q), 131 app_ctx->file_path, FILE_MODE_READ, 132 &viewer); CHKERRQ(ierr); 133 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 134 135 // Compute error with respect to reference solution 136 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 137 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 138 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 139 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 140 141 // Check error 142 if (error > app_ctx->test_tol) { 143 ierr = PetscPrintf(PETSC_COMM_WORLD, 144 "Test failed with error norm %g\n", 145 (double)error); CHKERRQ(ierr); 146 } 147 148 // Cleanup 149 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 150 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 151 152 PetscFunctionReturn(0); 153 } 154 155 // Get error for problems with exact solutions 156 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, 157 PetscScalar final_time) { 158 PetscInt loc_nodes; 159 Vec Q_exact, Q_exact_loc; 160 PetscReal rel_error, norm_error, norm_exact; 161 PetscErrorCode ierr; 162 PetscFunctionBegin; 163 164 // Get exact solution at final time 165 ierr = DMCreateGlobalVector(dm, &Q_exact); CHKERRQ(ierr); 166 ierr = DMGetLocalVector(dm, &Q_exact_loc); CHKERRQ(ierr); 167 ierr = VecGetSize(Q_exact_loc, &loc_nodes); CHKERRQ(ierr); 168 ierr = ICs_FixMultiplicity(dm, ceed_data, user, Q_exact_loc, Q_exact, 169 final_time); 170 CHKERRQ(ierr); 171 172 // Get |exact solution - obtained solution| 173 ierr = VecNorm(Q_exact, NORM_1, &norm_exact); CHKERRQ(ierr); 174 ierr = VecAXPY(Q, -1.0, Q_exact); CHKERRQ(ierr); 175 ierr = VecNorm(Q, NORM_1, &norm_error); CHKERRQ(ierr); 176 177 // Compute relative error 178 rel_error = norm_error / norm_exact; 179 180 // Output relative error 181 ierr = PetscPrintf(PETSC_COMM_WORLD, 182 "Relative Error: %g\n", 183 (double)rel_error); CHKERRQ(ierr); 184 // Cleanup 185 ierr = DMRestoreLocalVector(dm, &Q_exact_loc); CHKERRQ(ierr); 186 ierr = VecDestroy(&Q_exact); CHKERRQ(ierr); 187 188 PetscFunctionReturn(0); 189 } 190 191 // Post-processing 192 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, 193 ProblemData *problem, User user, 194 Vec Q, PetscScalar final_time) { 195 PetscInt steps; 196 PetscErrorCode ierr; 197 PetscFunctionBegin; 198 199 // Print relative error 200 if (problem->non_zero_time && !user->app_ctx->test_mode) { 201 ierr = GetError_NS(ceed_data, dm, user, Q, final_time); CHKERRQ(ierr); 202 } 203 204 // Print final time and number of steps 205 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 206 if (!user->app_ctx->test_mode) { 207 ierr = PetscPrintf(PETSC_COMM_WORLD, 208 "Time integrator took %" PetscInt_FMT " time steps to reach final time %g\n", 209 steps, (double)final_time); CHKERRQ(ierr); 210 } 211 212 // Output numerical values from command line 213 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 214 215 // Compare reference solution values with current test run for CI 216 if (user->app_ctx->test_mode) { 217 ierr = RegressionTests_NS(user->app_ctx, Q); CHKERRQ(ierr); 218 } 219 220 PetscFunctionReturn(0); 221 } 222 223 // Gather initial Q values in case of continuation of simulation 224 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q) { 225 226 PetscViewer viewer; 227 char file_path[PETSC_MAX_PATH_LEN]; 228 PetscErrorCode ierr; 229 PetscFunctionBegin; 230 231 // Read input 232 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 233 app_ctx->output_dir); CHKERRQ(ierr); 234 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 235 CHKERRQ(ierr); 236 237 // Load Q from existent solution 238 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 239 240 // Cleanup 241 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 242 243 PetscFunctionReturn(0); 244 } 245 246 // Record boundary values from initial condition 247 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc) { 248 249 Vec Qbc; 250 PetscErrorCode ierr; 251 PetscFunctionBegin; 252 253 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 254 ierr = VecCopy(Q_loc, Qbc); CHKERRQ(ierr); 255 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 256 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 257 ierr = VecAXPY(Qbc, -1., Q_loc); CHKERRQ(ierr); 258 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 259 ierr = PetscObjectComposeFunction((PetscObject)dm, 260 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 261 CHKERRQ(ierr); 262 263 PetscFunctionReturn(0); 264 } 265 266 // Free a plain data context that was allocated using PETSc; returning libCEED error codes 267 int FreeContextPetsc(void *data) { 268 if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, 269 "PetscFree failed"); 270 return CEED_ERROR_SUCCESS; 271 } 272