1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Miscellaneous utility functions 19 20 #include "../navierstokes.h" 21 22 PetscErrorCode ICs_FixMultiplicity(DM dm, CeedData ceed_data, Vec Q_loc, Vec Q, 23 CeedScalar time) { 24 PetscErrorCode ierr; 25 PetscFunctionBeginUser; 26 27 // --------------------------------------------------------------------------- 28 // Update SetupContext 29 // --------------------------------------------------------------------------- 30 SetupContext setup_ctx; 31 CeedQFunctionContextGetData(ceed_data->setup_context, CEED_MEM_HOST, 32 (void **)&setup_ctx); 33 setup_ctx->time = time; 34 CeedQFunctionContextRestoreData(ceed_data->setup_context, (void **)&setup_ctx); 35 36 // --------------------------------------------------------------------------- 37 // ICs 38 // --------------------------------------------------------------------------- 39 // -- CEED Restriction 40 CeedVector q0_ceed; 41 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q0_ceed, NULL); 42 43 // -- Place PETSc vector in CEED vector 44 CeedScalar *q0; 45 PetscMemType q0_mem_type; 46 ierr = VecGetArrayAndMemType(Q_loc, (PetscScalar **)&q0, &q0_mem_type); 47 CHKERRQ(ierr); 48 CeedVectorSetArray(q0_ceed, MemTypeP2C(q0_mem_type), CEED_USE_POINTER, q0); 49 50 // -- Apply CEED Operator 51 CeedOperatorApply(ceed_data->op_ics, ceed_data->x_coord, q0_ceed, 52 CEED_REQUEST_IMMEDIATE); 53 54 // -- Restore vectors 55 CeedVectorTakeArray(q0_ceed, MemTypeP2C(q0_mem_type), NULL); 56 ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q0); 57 CHKERRQ(ierr); 58 59 // -- Local-to-Global 60 ierr = VecZeroEntries(Q); CHKERRQ(ierr); 61 ierr = DMLocalToGlobal(dm, Q_loc, ADD_VALUES, Q); CHKERRQ(ierr); 62 63 // --------------------------------------------------------------------------- 64 // Fix multiplicity for output of ICs 65 // --------------------------------------------------------------------------- 66 // -- CEED Restriction 67 CeedVector mult_vec; 68 CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &mult_vec, NULL); 69 70 // -- Place PETSc vector in CEED vector 71 CeedScalar *mult; 72 PetscMemType m_mem_type; 73 Vec multiplicity_loc; 74 ierr = DMGetLocalVector(dm, &multiplicity_loc); CHKERRQ(ierr); 75 ierr = VecGetArrayAndMemType(multiplicity_loc, (PetscScalar **)&mult, 76 &m_mem_type); 77 CHKERRQ(ierr); 78 CeedVectorSetArray(mult_vec, MemTypeP2C(m_mem_type), CEED_USE_POINTER, mult); 79 CHKERRQ(ierr); 80 81 // -- Get multiplicity 82 CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, mult_vec); 83 84 // -- Restore vectors 85 CeedVectorTakeArray(mult_vec, MemTypeP2C(m_mem_type), NULL); 86 ierr = VecRestoreArrayReadAndMemType(multiplicity_loc, 87 (const PetscScalar **)&mult); CHKERRQ(ierr); 88 89 // -- Local-to-Global 90 Vec multiplicity; 91 ierr = DMGetGlobalVector(dm, &multiplicity); CHKERRQ(ierr); 92 ierr = VecZeroEntries(multiplicity); CHKERRQ(ierr); 93 ierr = DMLocalToGlobal(dm, multiplicity_loc, ADD_VALUES, multiplicity); 94 CHKERRQ(ierr); 95 96 // -- Fix multiplicity 97 ierr = VecPointwiseDivide(Q, Q, multiplicity); CHKERRQ(ierr); 98 ierr = VecPointwiseDivide(Q_loc, Q_loc, multiplicity_loc); CHKERRQ(ierr); 99 100 // -- Restore vectors 101 ierr = DMRestoreLocalVector(dm, &multiplicity_loc); CHKERRQ(ierr); 102 ierr = DMRestoreGlobalVector(dm, &multiplicity); CHKERRQ(ierr); 103 104 // Cleanup 105 CeedVectorDestroy(&mult_vec); 106 CeedVectorDestroy(&q0_ceed); 107 108 PetscFunctionReturn(0); 109 } 110 111 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, 112 PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 113 Vec cell_geom_FVM, Vec grad_FVM) { 114 115 Vec Qbc; 116 PetscErrorCode ierr; 117 PetscFunctionBegin; 118 119 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 120 ierr = VecAXPY(Q_loc, 1., Qbc); CHKERRQ(ierr); 121 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 122 123 PetscFunctionReturn(0); 124 } 125 126 // Compare reference solution values with current test run for CI 127 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) { 128 129 Vec Qref; 130 PetscViewer viewer; 131 PetscReal error, Qrefnorm; 132 PetscErrorCode ierr; 133 PetscFunctionBegin; 134 135 // Read reference file 136 ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr); 137 ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)Q), 138 app_ctx->file_path, FILE_MODE_READ, 139 &viewer); CHKERRQ(ierr); 140 ierr = VecLoad(Qref, viewer); CHKERRQ(ierr); 141 142 // Compute error with respect to reference solution 143 ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr); 144 ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr); 145 ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr); 146 ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr); 147 148 // Check error 149 if (error > app_ctx->test_tol) { 150 ierr = PetscPrintf(PETSC_COMM_WORLD, 151 "Test failed with error norm %g\n", 152 (double)error); CHKERRQ(ierr); 153 } 154 155 // Cleanup 156 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 157 ierr = VecDestroy(&Qref); CHKERRQ(ierr); 158 159 PetscFunctionReturn(0); 160 } 161 162 // Get error for problems with exact solutions 163 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, AppCtx app_ctx, Vec Q, 164 PetscScalar final_time) { 165 PetscInt loc_nodes; 166 Vec Q_exact, Q_exact_loc; 167 PetscReal rel_error, norm_error, norm_exact; 168 PetscErrorCode ierr; 169 PetscFunctionBegin; 170 171 // Get exact solution at final time 172 ierr = DMCreateGlobalVector(dm, &Q_exact); CHKERRQ(ierr); 173 ierr = DMGetLocalVector(dm, &Q_exact_loc); CHKERRQ(ierr); 174 ierr = VecGetSize(Q_exact_loc, &loc_nodes); CHKERRQ(ierr); 175 ierr = ICs_FixMultiplicity(dm, ceed_data, Q_exact_loc, Q_exact, final_time); 176 CHKERRQ(ierr); 177 178 // Get |exact solution - obtained solution| 179 ierr = VecNorm(Q_exact, NORM_1, &norm_exact); CHKERRQ(ierr); 180 ierr = VecAXPY(Q, -1.0, Q_exact); CHKERRQ(ierr); 181 ierr = VecNorm(Q, NORM_1, &norm_error); CHKERRQ(ierr); 182 183 // Compute relative error 184 rel_error = norm_error / norm_exact; 185 186 // Output relative error 187 ierr = PetscPrintf(PETSC_COMM_WORLD, 188 "Relative Error: %g\n", 189 (double)rel_error); CHKERRQ(ierr); 190 // Cleanup 191 ierr = DMRestoreLocalVector(dm, &Q_exact_loc); CHKERRQ(ierr); 192 ierr = VecDestroy(&Q_exact); CHKERRQ(ierr); 193 194 PetscFunctionReturn(0); 195 } 196 197 // Post-processing 198 PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, 199 ProblemData *problem, AppCtx app_ctx, 200 Vec Q, PetscScalar final_time) { 201 PetscInt steps; 202 PetscErrorCode ierr; 203 PetscFunctionBegin; 204 205 // Print relative error 206 if (problem->non_zero_time && !app_ctx->test_mode) { 207 ierr = GetError_NS(ceed_data, dm, app_ctx, Q, final_time); CHKERRQ(ierr); 208 } 209 210 // Print final time and number of steps 211 ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr); 212 if (!app_ctx->test_mode) { 213 ierr = PetscPrintf(PETSC_COMM_WORLD, 214 "Time integrator took %D time steps to reach final time %g\n", 215 steps, (double)final_time); CHKERRQ(ierr); 216 } 217 218 // Output numerical values from command line 219 ierr = VecViewFromOptions(Q, NULL, "-vec_view"); CHKERRQ(ierr); 220 221 // Compare reference solution values with current test run for CI 222 if (app_ctx->test_mode) { 223 ierr = RegressionTests_NS(app_ctx, Q); CHKERRQ(ierr); 224 } 225 226 PetscFunctionReturn(0); 227 } 228 229 // Gather initial Q values in case of continuation of simulation 230 PetscErrorCode SetupICsFromBinary(MPI_Comm comm, AppCtx app_ctx, Vec Q) { 231 232 PetscViewer viewer; 233 char file_path[PETSC_MAX_PATH_LEN]; 234 PetscErrorCode ierr; 235 PetscFunctionBegin; 236 237 // Read input 238 ierr = PetscSNPrintf(file_path, sizeof file_path, "%s/ns-solution.bin", 239 app_ctx->output_dir); CHKERRQ(ierr); 240 ierr = PetscViewerBinaryOpen(comm, file_path, FILE_MODE_READ, &viewer); 241 CHKERRQ(ierr); 242 243 // Load Q from existent solution 244 ierr = VecLoad(Q, viewer); CHKERRQ(ierr); 245 246 // Cleanup 247 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 248 249 PetscFunctionReturn(0); 250 } 251 252 // Record boundary values from initial condition 253 PetscErrorCode SetBCsFromICs_NS(DM dm, Vec Q, Vec Q_loc) { 254 255 Vec Qbc; 256 PetscErrorCode ierr; 257 PetscFunctionBegin; 258 259 ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 260 ierr = VecCopy(Q_loc, Qbc); CHKERRQ(ierr); 261 ierr = VecZeroEntries(Q_loc); CHKERRQ(ierr); 262 ierr = DMGlobalToLocal(dm, Q, INSERT_VALUES, Q_loc); CHKERRQ(ierr); 263 ierr = VecAXPY(Qbc, -1., Q_loc); CHKERRQ(ierr); 264 ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr); 265 ierr = PetscObjectComposeFunction((PetscObject)dm, 266 "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_NS); 267 CHKERRQ(ierr); 268 269 PetscFunctionReturn(0); 270 } 271