Line data Source code
1 : /********************************************************************/
2 : /* SOFTWARE COPYRIGHT NOTIFICATION */
3 : /* Cardinal */
4 : /* */
5 : /* (c) 2021 UChicago Argonne, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by UChicago Argonne, LLC */
9 : /* Under Contract No. DE-AC02-06CH11357 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* Prepared by Battelle Energy Alliance, LLC */
13 : /* Under Contract No. DE-AC07-05ID14517 */
14 : /* With the U. S. Department of Energy */
15 : /* */
16 : /* See LICENSE for full restrictions */
17 : /********************************************************************/
18 :
19 : #ifdef ENABLE_NEK_COUPLING
20 :
21 : #include "NekBoundaryFlux.h"
22 :
23 : registerMooseObject("CardinalApp", NekBoundaryFlux);
24 :
25 : extern nekrs::usrwrkIndices indices;
26 :
27 : InputParameters
28 492 : NekBoundaryFlux::validParams()
29 : {
30 492 : auto params = ConservativeFieldTransfer::validParams();
31 984 : params.addParam<Real>(
32 : "initial_flux_integral",
33 984 : 0,
34 : "Initial value to use for the 'postprocessor_to_conserve'; this initial value will be "
35 : "overridden once the coupled app executes its transfer of the boundary flux term integral "
36 : "into the 'postprocessor_to_conserve'. You may want to use this parameter if NekRS runs "
37 : "first, or if you are running NekRS in isolation but still want to apply a boundary flux "
38 : "term via Cardinal. Remember that this parameter is only used to normalize the flux, so you "
39 : "will need to populate an initial shape if you want to see this parameter take effect.");
40 :
41 984 : params.addParam<bool>(
42 : "conserve_flux_by_sideset",
43 984 : false,
44 : "Whether to conserve the flux by individual sideset (as opposed to lumping all sidesets "
45 : "together). Setting this option to true requires syntax changes in the input file to use "
46 : "vector postprocessors, and places restrictions on how the sidesets are set up.");
47 :
48 492 : params.addClassDescription("Reads/writes boundary flux data between NekRS and MOOSE.");
49 492 : return params;
50 0 : }
51 :
52 248 : NekBoundaryFlux::NekBoundaryFlux(const InputParameters & parameters)
53 : : ConservativeFieldTransfer(parameters),
54 246 : _conserve_flux_by_sideset(getParam<bool>("conserve_flux_by_sideset")),
55 492 : _initial_flux_integral(getParam<Real>("initial_flux_integral")),
56 246 : _boundary(_nek_mesh->boundary()),
57 494 : _reference_flux_integral(nekrs::referenceArea() * nekrs::nondimensionalDivisor(field::flux))
58 : {
59 246 : if (_direction == "to_nek")
60 : {
61 246 : if (_usrwrk_slot.size() > 1)
62 2 : paramError("usrwrk_slot",
63 : "'usrwrk_slot' must be of length 1 for boundary flux transfers; you have entered "
64 1 : "a vector of length " +
65 0 : Moose::stringify(_usrwrk_slot.size()));
66 :
67 489 : addExternalVariable(_usrwrk_slot[0], _variable);
68 244 : indices.flux = _usrwrk_slot[0] * nekrs::fieldOffset();
69 : }
70 :
71 244 : if (!_boundary)
72 1 : mooseError("NekBoundaryFlux can only be used when there is boundary coupling of NekRS with "
73 : "MOOSE, i.e. when 'boundary' is provided in NekRSMesh.");
74 :
75 243 : if (_direction == "from_nek")
76 0 : paramError("direction",
77 : "The NekBoundaryFlux currently only supports transfers 'to_nek'; contact the "
78 : "Cardinal developer team if you require reading of NekRS boundary flux terms.");
79 :
80 : // Check that the correct flux boundary condition is set on all of nekRS's
81 : // boundaries. To avoid throwing this error for test cases where we have a
82 : // [TEMPERATURE] block but set its solve to 'none', we screen by whether the
83 : // temperature solve is enabled
84 243 : if (_boundary && nekrs::hasTemperatureSolve())
85 437 : for (const auto & b : *_boundary)
86 219 : if (!nekrs::isHeatFluxBoundary(b))
87 1 : mooseError("In order to send a boundary heat flux to NekRS, you must have a heat flux "
88 1 : "condition for each 'boundary' set in 'NekRSMesh'! Boundary " +
89 2 : std::to_string(b) + " is of type '" + nekrs::temperatureBoundaryType(b) +
90 : "' instead of 'fixedGradient'.");
91 :
92 242 : if (!nekrs::hasTemperatureVariable())
93 1 : mooseError("In order to send a boundary heat flux to NekRS, your case files must have a "
94 1 : "[TEMPERATURE] block. Note that you can set 'solver = none' in '" +
95 1 : _nek_problem.casename() + ".par' if you don't want to solve for temperature.");
96 :
97 241 : if (!nekrs::hasTemperatureSolve())
98 68 : mooseWarning("By setting 'solver = none' for temperature in '" + _nek_problem.casename() +
99 : ".par', NekRS will not solve for temperature. The heat flux sent by this object "
100 : "will be unused.");
101 :
102 : // add the postprocessor that receives the flux integral for normalization
103 240 : if (_conserve_flux_by_sideset)
104 : {
105 40 : if (isParamSetByUser("initial_flux_integral"))
106 0 : mooseWarning("The 'initial_flux_integral' capability is not yet supported when "
107 : "'conserve_flux_by_sideset' is enabled. Please contact a Cardinal developer "
108 : "if this is hindering your use case.");
109 :
110 20 : auto vpp_params = _factory.getValidParams("ConstantVectorPostprocessor");
111 :
112 : // create zero initial values
113 20 : std::vector<std::vector<Real>> dummy_vals(1, std::vector<Real>(_boundary->size()));
114 20 : vpp_params.set<std::vector<std::vector<Real>>>("value") = dummy_vals;
115 20 : _nek_problem.addVectorPostprocessor(
116 20 : "ConstantVectorPostprocessor", _postprocessor_name, vpp_params);
117 20 : }
118 : else
119 440 : addExternalPostprocessor(_postprocessor_name, _initial_flux_integral);
120 :
121 240 : if (_conserve_flux_by_sideset)
122 20 : _flux_integral_vpp =
123 40 : &_nek_problem.getVectorPostprocessorValueByName(_postprocessor_name, "value");
124 : else
125 220 : _flux_integral = &getPostprocessorValueByName(_postprocessor_name);
126 :
127 240 : _flux_face = (double *)calloc(_n_per_surf, sizeof(double));
128 240 : _flux_elem = (double *)calloc(_n_per_vol, sizeof(double));
129 240 : }
130 :
131 464 : NekBoundaryFlux::~NekBoundaryFlux()
132 : {
133 232 : freePointer(_flux_face);
134 232 : freePointer(_flux_elem);
135 464 : }
136 :
137 : void
138 13314 : NekBoundaryFlux::sendDataToNek()
139 : {
140 26628 : _console << "Sending flux to NekRS boundary " << Moose::stringify(*_boundary) << "..."
141 13314 : << std::endl;
142 13314 : auto d = nekrs::nondimensionalDivisor(field::flux);
143 13314 : auto a = nekrs::nondimensionalAdditive(field::flux);
144 :
145 13314 : if (!_nek_mesh->volume())
146 : {
147 1863100 : for (unsigned int e = 0; e < _nek_mesh->numSurfaceElems(); e++)
148 : {
149 : // We can only write into the nekRS scratch space if that face is "owned" by the current
150 : // process
151 1850608 : if (nekrs::commRank() != _nek_mesh->boundaryCoupling().processor_id(e))
152 1561840 : continue;
153 :
154 288768 : _nek_problem.mapFaceDataToNekFace(e, _variable_number[_variable], d, a, &_flux_face);
155 288768 : _nek_problem.writeBoundarySolution(e, field::flux, _flux_face);
156 : }
157 : }
158 : else
159 : {
160 2975814 : for (unsigned int e = 0; e < _nek_mesh->numVolumeElems(); ++e)
161 : {
162 : // We can only write into the nekRS scratch space if that face is "owned" by the current
163 : // process
164 2974992 : if (nekrs::commRank() != _nek_mesh->volumeCoupling().processor_id(e))
165 2338656 : continue;
166 :
167 636336 : _nek_problem.mapFaceDataToNekVolume(e, _variable_number[_variable], d, a, &_flux_elem);
168 636336 : _nek_problem.writeVolumeSolution(e, field::flux, _flux_elem);
169 : }
170 : }
171 :
172 : // Because the NekRSMesh may be quite different from that used in the app solving for
173 : // the heat flux, we will need to normalize the flux on the nekRS side by the
174 : // flux computed by the coupled MOOSE app. For this and the next check of the
175 : // flux integral, we need to scale the integral back up again to the dimensional form
176 : // for the sake of comparison.
177 13314 : const Real scale_squared = _nek_mesh->scaling() * _nek_mesh->scaling();
178 13314 : const double nek_flux_print_mult = scale_squared * nekrs::nondimensionalDivisor(field::flux);
179 :
180 : // integrate the flux over each individual boundary
181 : std::vector<double> nek_flux_sidesets =
182 13314 : nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
183 :
184 : bool successful_normalization;
185 13314 : double normalized_nek_flux = 0.0;
186 :
187 : double total_moose_flux;
188 :
189 13314 : if (_conserve_flux_by_sideset)
190 : {
191 254 : auto moose_flux = *_flux_integral_vpp;
192 254 : if (moose_flux.size() != _boundary->size())
193 1 : mooseError("The sideset flux reporter transferred to NekRS must have a length equal to the "
194 : "number of entries in 'boundary'! Please check the values written to the "
195 : "'flux_integral' vector postprocessor.\n\n"
196 : "Length of reporter: ",
197 1 : moose_flux.size(),
198 : "\n",
199 : "Length of 'boundary': ",
200 1 : _boundary->size());
201 :
202 549 : for (std::size_t b = 0; b < _boundary->size(); ++b)
203 : {
204 592 : _console << "[boundary " << Moose::stringify((*_boundary)[b])
205 296 : << "]: Normalizing NekRS flux of "
206 592 : << Moose::stringify(nek_flux_sidesets[b] * nek_flux_print_mult)
207 592 : << " to the conserved MOOSE value of " << Moose::stringify(moose_flux[b])
208 296 : << std::endl;
209 :
210 296 : checkInitialFluxValues(nek_flux_sidesets[b], moose_flux[b]);
211 : }
212 :
213 253 : total_moose_flux = std::accumulate(moose_flux.begin(), moose_flux.end(), 0.0);
214 :
215 : // For the sake of printing diagnostics to the screen regarding the flux normalization,
216 : // we first scale the nek flux by any unit changes and then by the reference flux.
217 : successful_normalization =
218 253 : normalizeFluxBySideset(moose_flux, nek_flux_sidesets, normalized_nek_flux);
219 : }
220 : else
221 : {
222 13060 : auto moose_flux = *_flux_integral;
223 : const double nek_flux =
224 13060 : std::accumulate(nek_flux_sidesets.begin(), nek_flux_sidesets.end(), 0.0);
225 :
226 26120 : _console << "[boundary " << Moose::stringify(*_boundary)
227 13060 : << "]: Normalizing total NekRS flux of "
228 26120 : << Moose::stringify(nek_flux * nek_flux_print_mult)
229 13060 : << " to the conserved MOOSE value of " << Moose::stringify(moose_flux) << std::endl;
230 13060 : _console << Moose::stringify(nek_flux) << " " << nek_flux_print_mult << std::endl;
231 :
232 13060 : checkInitialFluxValues(nek_flux, moose_flux);
233 :
234 13060 : total_moose_flux = moose_flux;
235 :
236 : // For the sake of printing diagnostics to the screen regarding the flux normalization,
237 13060 : successful_normalization = normalizeFlux(moose_flux, nek_flux, normalized_nek_flux);
238 : }
239 :
240 13313 : if (!successful_normalization)
241 1 : mooseError(
242 : "Flux normalization process failed! NekRS integrated flux: ",
243 : normalized_nek_flux,
244 : " MOOSE integrated flux: ",
245 : total_moose_flux,
246 1 : normalizationHint(),
247 : "\n\n- If you set 'conserve_flux_by_sideset = true' and nodes are SHARED by boundaries "
248 : "(like on corners between sidesets), you will end up renormalizing those shared nodes "
249 : "once per sideset that they lie on. There is no guarantee that the total imposed flux "
250 : "would be preserved.");
251 13312 : }
252 :
253 : void
254 13356 : NekBoundaryFlux::checkInitialFluxValues(const Real & nek_flux, const Real & moose_flux) const
255 : {
256 13356 : const Real scale_squared = _nek_mesh->scaling() * _nek_mesh->scaling();
257 13356 : const double nek_flux_print_mult = scale_squared * nekrs::nondimensionalDivisor(field::flux);
258 :
259 : // If before normalization, there is a large difference between the nekRS imposed flux
260 : // and the MOOSE flux, this could mean that there is a poor match between the domains,
261 : // even if neither value is zero. For instance, if you forgot that the nekRS mesh is in
262 : // units of centimeters, but you're coupling to an app based in meters, the fluxes will
263 : // be very different from one another.
264 13356 : if (moose_flux && (std::abs(nek_flux * nek_flux_print_mult - moose_flux) / moose_flux) > 0.25)
265 72 : mooseDoOnce(mooseWarning(
266 : "NekRS flux differs from MOOSE flux by more than 25\%! This is NOT necessarily a problem - "
267 : "but it could indicate that your geometries don't line up properly or something is amiss "
268 : "with your transfer. We recommend opening the output files to visually inspect the flux in "
269 : "both the main and sub applications to check that the fields look correct."));
270 13356 : }
271 :
272 : bool
273 253 : NekBoundaryFlux::normalizeFluxBySideset(const std::vector<double> & moose_integral,
274 : std::vector<double> & nek_integral,
275 : double & normalized_nek_integral)
276 : {
277 : // scale the nek flux to dimensional form for the sake of normalizing against
278 : // a dimensional MOOSE flux
279 549 : for (auto & i : nek_integral)
280 296 : i *= _reference_flux_integral;
281 :
282 253 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
283 253 : mesh_t * mesh = nekrs::temperatureMesh();
284 253 : auto nek_boundary_coupling = _nek_mesh->boundaryCoupling();
285 :
286 146405 : for (int k = 0; k < nek_boundary_coupling.total_n_faces; ++k)
287 : {
288 146152 : if (nek_boundary_coupling.process[k] == nekrs::commRank())
289 : {
290 25264 : int i = nek_boundary_coupling.element[k];
291 25264 : int j = nek_boundary_coupling.face[k];
292 :
293 25264 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
294 25264 : auto it = std::find(_boundary->begin(), _boundary->end(), face_id);
295 : auto b_index = it - _boundary->begin();
296 :
297 : // avoid divide-by-zero
298 : double ratio = 1.0;
299 25264 : if (std::abs(nek_integral[b_index]) > _abs_tol)
300 24512 : ratio = moose_integral[b_index] / nek_integral[b_index];
301 :
302 25264 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
303 :
304 670768 : for (int v = 0; v < mesh->Nfp; ++v)
305 : {
306 645504 : int id = mesh->vmapM[offset + v];
307 645504 : nrs->usrwrk[indices.flux + id] *= ratio;
308 : }
309 : }
310 : }
311 :
312 : // check that the normalization worked properly - confirm against dimensional form
313 253 : auto integrals = nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
314 253 : normalized_nek_integral =
315 253 : std::accumulate(integrals.begin(), integrals.end(), 0.0) * _reference_flux_integral;
316 : double total_moose_integral = std::accumulate(moose_integral.begin(), moose_integral.end(), 0.0);
317 : bool low_rel_err =
318 253 : std::abs(total_moose_integral) > _abs_tol
319 253 : ? std::abs(normalized_nek_integral - total_moose_integral) / total_moose_integral <
320 249 : _rel_tol
321 : : true;
322 253 : bool low_abs_err = std::abs(normalized_nek_integral - total_moose_integral) < _abs_tol;
323 :
324 506 : return low_rel_err || low_abs_err;
325 253 : }
326 :
327 : bool
328 13060 : NekBoundaryFlux::normalizeFlux(const double moose_integral,
329 : double nek_integral,
330 : double & normalized_nek_integral)
331 : {
332 : // scale the nek flux to dimensional form for the sake of normalizing against
333 : // a dimensional MOOSE flux
334 13060 : nek_integral *= _reference_flux_integral;
335 :
336 13060 : auto nek_boundary_coupling = _nek_mesh->boundaryCoupling();
337 :
338 : // avoid divide-by-zero
339 13060 : if (std::abs(nek_integral) < _abs_tol)
340 : return true;
341 :
342 12876 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
343 12876 : mesh_t * mesh = nekrs::temperatureMesh();
344 :
345 12876 : const double ratio = moose_integral / nek_integral;
346 :
347 2049276 : for (int k = 0; k < nek_boundary_coupling.total_n_faces; ++k)
348 : {
349 2036400 : if (nek_boundary_coupling.process[k] == nekrs::commRank())
350 : {
351 390624 : int i = nek_boundary_coupling.element[k];
352 390624 : int j = nek_boundary_coupling.face[k];
353 390624 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
354 :
355 7235064 : for (int v = 0; v < mesh->Nfp; ++v)
356 : {
357 6844440 : int id = mesh->vmapM[offset + v];
358 6844440 : nrs->usrwrk[indices.flux + id] *= ratio;
359 : }
360 : }
361 : }
362 :
363 : // check that the normalization worked properly - confirm against dimensional form
364 12876 : auto integrals = nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
365 12876 : normalized_nek_integral =
366 12876 : std::accumulate(integrals.begin(), integrals.end(), 0.0) * _reference_flux_integral;
367 12876 : bool low_rel_err = std::abs(normalized_nek_integral - moose_integral) / moose_integral < _rel_tol;
368 12876 : bool low_abs_err = std::abs(normalized_nek_integral - moose_integral) < _abs_tol;
369 :
370 12876 : return low_rel_err || low_abs_err;
371 13060 : }
372 :
373 : #endif
|
