xref: /libCEED/examples/fluids/src/misc.c (revision 88626eed6564cd43033d3137230605fb5f962840) 
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, Vec Q_loc, Vec Q,
14                                    CeedScalar time) {
15   PetscErrorCode ierr;
16   PetscFunctionBeginUser;
17 
18   // ---------------------------------------------------------------------------
19   // Update SetupContext
20   // ---------------------------------------------------------------------------
21   SetupContext setup_ctx;
22   CeedQFunctionContextGetData(ceed_data->setup_context, CEED_MEM_HOST,
23                               (void **)&setup_ctx);
24   setup_ctx->time = time;
25   CeedQFunctionContextRestoreData(ceed_data->setup_context, (void **)&setup_ctx);
26 
27   // ---------------------------------------------------------------------------
28   // ICs
29   // ---------------------------------------------------------------------------
30   // -- CEED Restriction
31   CeedVector q0_ceed;
32   CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q0_ceed, NULL);
33 
34   // -- Place PETSc vector in CEED vector
35   CeedScalar *q0;
36   PetscMemType q0_mem_type;
37   ierr = VecGetArrayAndMemType(Q_loc, (PetscScalar **)&q0, &q0_mem_type);
38   CHKERRQ(ierr);
39   CeedVectorSetArray(q0_ceed, MemTypeP2C(q0_mem_type), CEED_USE_POINTER, q0);
40 
41   // -- Apply CEED Operator
42   CeedOperatorApply(ceed_data->op_ics, ceed_data->x_coord, q0_ceed,
43                     CEED_REQUEST_IMMEDIATE);
44 
45   // -- Restore vectors
46   CeedVectorTakeArray(q0_ceed, MemTypeP2C(q0_mem_type), NULL);
47   ierr = VecRestoreArrayReadAndMemType(Q_loc, (const PetscScalar **)&q0);
48   CHKERRQ(ierr);
49 
50   // -- Local-to-Global
51   ierr = VecZeroEntries(Q); CHKERRQ(ierr);
52   ierr = DMLocalToGlobal(dm, Q_loc, ADD_VALUES, Q); CHKERRQ(ierr);
53 
54   // ---------------------------------------------------------------------------
55   // Fix multiplicity for output of ICs
56   // ---------------------------------------------------------------------------
57   // -- CEED Restriction
58   CeedVector mult_vec;
59   CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &mult_vec, NULL);
60 
61   // -- Place PETSc vector in CEED vector
62   CeedScalar *mult;
63   PetscMemType m_mem_type;
64   Vec multiplicity_loc;
65   ierr = DMGetLocalVector(dm, &multiplicity_loc); CHKERRQ(ierr);
66   ierr = VecGetArrayAndMemType(multiplicity_loc, (PetscScalar **)&mult,
67                                &m_mem_type);
68   CHKERRQ(ierr);
69   CeedVectorSetArray(mult_vec, MemTypeP2C(m_mem_type), CEED_USE_POINTER, mult);
70   CHKERRQ(ierr);
71 
72   // -- Get multiplicity
73   CeedElemRestrictionGetMultiplicity(ceed_data->elem_restr_q, mult_vec);
74 
75   // -- Restore vectors
76   CeedVectorTakeArray(mult_vec, MemTypeP2C(m_mem_type), NULL);
77   ierr = VecRestoreArrayReadAndMemType(multiplicity_loc,
78                                        (const PetscScalar **)&mult); CHKERRQ(ierr);
79 
80   // -- Local-to-Global
81   Vec multiplicity;
82   ierr = DMGetGlobalVector(dm, &multiplicity); CHKERRQ(ierr);
83   ierr = VecZeroEntries(multiplicity); CHKERRQ(ierr);
84   ierr = DMLocalToGlobal(dm, multiplicity_loc, ADD_VALUES, multiplicity);
85   CHKERRQ(ierr);
86 
87   // -- Fix multiplicity
88   ierr = VecPointwiseDivide(Q, Q, multiplicity); CHKERRQ(ierr);
89   ierr = VecPointwiseDivide(Q_loc, Q_loc, multiplicity_loc); CHKERRQ(ierr);
90 
91   // -- Restore vectors
92   ierr = DMRestoreLocalVector(dm, &multiplicity_loc); CHKERRQ(ierr);
93   ierr = DMRestoreGlobalVector(dm, &multiplicity); CHKERRQ(ierr);
94 
95   // Cleanup
96   CeedVectorDestroy(&mult_vec);
97   CeedVectorDestroy(&q0_ceed);
98 
99   PetscFunctionReturn(0);
100 }
101 
102 PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm,
103     PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM,
104     Vec cell_geom_FVM, Vec grad_FVM) {
105 
106   Vec            Qbc;
107   PetscErrorCode ierr;
108   PetscFunctionBegin;
109 
110   ierr = DMGetNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr);
111   ierr = VecAXPY(Q_loc, 1., Qbc); CHKERRQ(ierr);
112   ierr = DMRestoreNamedLocalVector(dm, "Qbc", &Qbc); CHKERRQ(ierr);
113 
114   PetscFunctionReturn(0);
115 }
116 
117 // Compare reference solution values with current test run for CI
118 PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) {
119 
120   Vec            Qref;
121   PetscViewer    viewer;
122   PetscReal      error, Qrefnorm;
123   PetscErrorCode ierr;
124   PetscFunctionBegin;
125 
126   // Read reference file
127   ierr = VecDuplicate(Q, &Qref); CHKERRQ(ierr);
128   ierr = PetscViewerBinaryOpen(PetscObjectComm((PetscObject)Q),
129                                app_ctx->file_path, FILE_MODE_READ,
130                                &viewer); CHKERRQ(ierr);
131   ierr = VecLoad(Qref, viewer); CHKERRQ(ierr);
132 
133   // Compute error with respect to reference solution
134   ierr = VecAXPY(Q, -1.0, Qref); CHKERRQ(ierr);
135   ierr = VecNorm(Qref, NORM_MAX, &Qrefnorm); CHKERRQ(ierr);
136   ierr = VecScale(Q, 1./Qrefnorm); CHKERRQ(ierr);
137   ierr = VecNorm(Q, NORM_MAX, &error); CHKERRQ(ierr);
138 
139   // Check error
140   if (error > app_ctx->test_tol) {
141     ierr = PetscPrintf(PETSC_COMM_WORLD,
142                        "Test failed with error norm %g\n",
143                        (double)error); CHKERRQ(ierr);
144   }
145 
146   // Cleanup
147   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
148   ierr = VecDestroy(&Qref); CHKERRQ(ierr);
149 
150   PetscFunctionReturn(0);
151 }
152 
153 // Get error for problems with exact solutions
154 PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, AppCtx app_ctx, Vec Q,
155                            PetscScalar final_time) {
156   PetscInt       loc_nodes;
157   Vec            Q_exact, Q_exact_loc;
158   PetscReal      rel_error, norm_error, norm_exact;
159   PetscErrorCode ierr;
160   PetscFunctionBegin;
161 
162   // Get exact solution at final time
163   ierr = DMCreateGlobalVector(dm, &Q_exact); CHKERRQ(ierr);
164   ierr = DMGetLocalVector(dm, &Q_exact_loc); CHKERRQ(ierr);
165   ierr = VecGetSize(Q_exact_loc, &loc_nodes); CHKERRQ(ierr);
166   ierr = ICs_FixMultiplicity(dm, ceed_data, Q_exact_loc, Q_exact, 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, AppCtx app_ctx,
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 && !app_ctx->test_mode) {
198     ierr = GetError_NS(ceed_data, dm, app_ctx, Q, final_time); CHKERRQ(ierr);
199   }
200 
201   // Print final time and number of steps
202   ierr = TSGetStepNumber(ts, &steps); CHKERRQ(ierr);
203   if (!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 (app_ctx->test_mode) {
214     ierr = RegressionTests_NS(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