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