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 525 : NekBoundaryFlux::validParams()
29 : {
30 525 : auto params = ConservativeFieldTransfer::validParams();
31 1050 : params.addParam<Real>(
32 : "initial_flux_integral",
33 1050 : 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 1050 : params.addParam<bool>(
42 : "conserve_flux_by_sideset",
43 1050 : 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 525 : params.addClassDescription("Reads/writes boundary flux data between NekRS and MOOSE.");
49 525 : return params;
50 0 : }
51 :
52 265 : NekBoundaryFlux::NekBoundaryFlux(const InputParameters & parameters)
53 : : ConservativeFieldTransfer(parameters),
54 263 : _conserve_flux_by_sideset(getParam<bool>("conserve_flux_by_sideset")),
55 526 : _initial_flux_integral(getParam<Real>("initial_flux_integral")),
56 263 : _boundary(_nek_mesh->boundary()),
57 528 : _reference_flux_integral(nekrs::referenceArea() * nekrs::nondimensionalDivisor(field::flux))
58 : {
59 263 : if (!_boundary)
60 1 : mooseError("NekBoundaryFlux can only be used when there is boundary coupling of NekRS with "
61 : "MOOSE, i.e. when 'boundary' is provided in NekRSMesh.");
62 :
63 262 : if (!nekrs::hasTemperatureVariable())
64 1 : mooseError("In order to read or write NekRS's boundary heat flux, your case files must have a "
65 1 : "[TEMPERATURE] block. Note that you can set 'solver = none' in '" +
66 1 : _nek_problem.casename() + ".par' if you don't want to solve for temperature.");
67 :
68 : // add the variables for the coupling and perform checks on problem setup
69 261 : if (_direction == "from_nek")
70 : {
71 17 : if (_conserve_flux_by_sideset)
72 0 : paramError("conserve_flux_by_sideset",
73 : "When 'direction = from_nek', the 'conserve_flux_by_sideset' option is not yet "
74 : "supported. Contact the Cardinal developer team if you need this feature.");
75 :
76 34 : checkUnusedParam(parameters, "initial_flux_integral", "'direction = from_nek'");
77 17 : addExternalVariable(_variable);
78 :
79 : // right now, all of our systems used for transferring data assume that if we have a volume
80 : // mesh mirror, that we will be writing data to that entire mesh mirror. But this is not the
81 : // case if we want to write a heat flux - since the notion of a unit outward normal is s
82 : // surface quantity. For now, just prevent users from doing this.
83 17 : if (_nek_mesh->volume())
84 1 : mooseError("The NekBoundaryFlux does not currently support writing heat flux on a boundary "
85 : "when the Mesh has 'volume = true.' Please contact the Cardinal developer team if "
86 : "you require this feature.");
87 : }
88 : else
89 : {
90 244 : if (_usrwrk_slot.size() > 1)
91 2 : paramError("usrwrk_slot",
92 : "'usrwrk_slot' must be of length 1 for boundary flux transfers; you have entered "
93 1 : "a vector of length " +
94 0 : Moose::stringify(_usrwrk_slot.size()));
95 :
96 485 : addExternalVariable(_usrwrk_slot[0], _variable);
97 242 : indices.flux = _usrwrk_slot[0] * nekrs::fieldOffset();
98 :
99 : // Check that the correct flux boundary condition is set on all of nekRS's
100 : // boundaries. To avoid throwing this error for test cases where we have a
101 : // [TEMPERATURE] block but set its solve to 'none', we screen by whether the
102 : // temperature solve is enabled
103 242 : if (_boundary && nekrs::hasTemperatureSolve())
104 437 : for (const auto & b : *_boundary)
105 219 : if (!nekrs::isHeatFluxBoundary(b))
106 1 : mooseError("In order to send a boundary heat flux to NekRS, you must have a heat flux "
107 1 : "condition for each 'boundary' set in 'NekRSMesh'! Boundary " +
108 2 : std::to_string(b) + " is of type '" + nekrs::temperatureBoundaryType(b) +
109 : "' instead of 'fixedGradient'.");
110 :
111 241 : if (!nekrs::hasTemperatureSolve())
112 68 : mooseWarning("By setting 'solver = none' for temperature in '" + _nek_problem.casename() +
113 : ".par', NekRS will not solve for temperature. The heat flux sent by this object "
114 : "will be unused.");
115 : }
116 :
117 : // add postprocessors to hold integrated heat flux
118 256 : if (_direction == "to_nek")
119 : {
120 : // add the postprocessor that receives the flux integral for normalization
121 240 : if (_conserve_flux_by_sideset)
122 : {
123 40 : if (isParamSetByUser("initial_flux_integral"))
124 0 : mooseWarning("The 'initial_flux_integral' capability is not yet supported when "
125 : "'conserve_flux_by_sideset' is enabled. Please contact a Cardinal developer "
126 : "if this is hindering your use case.");
127 :
128 20 : auto vpp_params = _factory.getValidParams("ConstantVectorPostprocessor");
129 :
130 : // create zero initial values
131 20 : std::vector<std::vector<Real>> dummy_vals(1, std::vector<Real>(_boundary->size()));
132 20 : vpp_params.set<std::vector<std::vector<Real>>>("value") = dummy_vals;
133 20 : _nek_problem.addVectorPostprocessor(
134 20 : "ConstantVectorPostprocessor", _postprocessor_name, vpp_params);
135 20 : }
136 : else
137 440 : addExternalPostprocessor(_postprocessor_name, _initial_flux_integral);
138 :
139 240 : if (_conserve_flux_by_sideset)
140 20 : _flux_integral_vpp =
141 40 : &_nek_problem.getVectorPostprocessorValueByName(_postprocessor_name, "value");
142 : else
143 220 : _flux_integral = &getPostprocessorValueByName(_postprocessor_name);
144 : }
145 : else
146 : {
147 : // create a NekHeatFluxIntegral postprocessor to compute the heat flux
148 16 : auto pp_params = _factory.getValidParams("NekHeatFluxIntegral");
149 16 : pp_params.set<std::vector<int>>("boundary") = *_boundary;
150 :
151 : // we do not need to check for duplicate names, because MOOSE already handles
152 : // this error checking
153 16 : _nek_problem.addPostprocessor("NekHeatFluxIntegral", _postprocessor_name, pp_params);
154 16 : _flux_integral = &getPostprocessorValueByName(_postprocessor_name);
155 16 : }
156 256 : }
157 :
158 : void
159 244 : NekBoundaryFlux::readDataFromNek()
160 : {
161 244 : if (!_nek_mesh->volume())
162 244 : _nek_problem.boundarySolution(field::flux, _external_data);
163 : else
164 0 : _nek_problem.volumeSolution(field::flux, _external_data);
165 :
166 244 : fillAuxVariable(_variable_number[_variable], _external_data);
167 :
168 : // the postprocessor is automatically populated because its execution controls its writing
169 244 : }
170 :
171 : void
172 13314 : NekBoundaryFlux::sendDataToNek()
173 : {
174 26628 : _console << "Sending flux to NekRS boundary " << Moose::stringify(*_boundary) << "..."
175 13314 : << std::endl;
176 13314 : auto d = nekrs::nondimensionalDivisor(field::flux);
177 13314 : auto a = nekrs::nondimensionalAdditive(field::flux);
178 :
179 13314 : if (!_nek_mesh->volume())
180 : {
181 1863100 : for (unsigned int e = 0; e < _nek_mesh->numSurfaceElems(); e++)
182 : {
183 : // We can only write into the nekRS scratch space if that face is "owned" by the current
184 : // process
185 1850608 : if (nekrs::commRank() != _nek_mesh->boundaryCoupling().processor_id(e))
186 1561840 : continue;
187 :
188 288768 : _nek_problem.mapFaceDataToNekFace(e, _variable_number[_variable], d, a, &_v_face);
189 288768 : _nek_problem.writeBoundarySolution(e, field::flux, _v_face);
190 : }
191 : }
192 : else
193 : {
194 2975814 : for (unsigned int e = 0; e < _nek_mesh->numVolumeElems(); ++e)
195 : {
196 : // We can only write into the nekRS scratch space if that face is "owned" by the current
197 : // process
198 2974992 : if (nekrs::commRank() != _nek_mesh->volumeCoupling().processor_id(e))
199 2338656 : continue;
200 :
201 636336 : _nek_problem.mapFaceDataToNekVolume(e, _variable_number[_variable], d, a, &_v_elem);
202 636336 : _nek_problem.writeVolumeSolution(e, field::flux, _v_elem);
203 : }
204 : }
205 :
206 : // Because the NekRSMesh may be quite different from that used in the app solving for
207 : // the heat flux, we will need to normalize the flux on the nekRS side by the
208 : // flux computed by the coupled MOOSE app. For this and the next check of the
209 : // flux integral, we need to scale the integral back up again to the dimensional form
210 : // for the sake of comparison.
211 13314 : const Real scale_squared = _nek_mesh->scaling() * _nek_mesh->scaling();
212 13314 : const double nek_flux_print_mult = scale_squared * nekrs::nondimensionalDivisor(field::flux);
213 :
214 : // integrate the flux over each individual boundary
215 : std::vector<double> nek_flux_sidesets =
216 13314 : nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
217 :
218 : bool successful_normalization;
219 13314 : double normalized_nek_flux = 0.0;
220 :
221 : double total_moose_flux;
222 :
223 13314 : if (_conserve_flux_by_sideset)
224 : {
225 254 : auto moose_flux = *_flux_integral_vpp;
226 254 : if (moose_flux.size() != _boundary->size())
227 1 : mooseError("The sideset flux reporter transferred to NekRS must have a length equal to the "
228 : "number of entries in 'boundary'! Please check the values written to the "
229 : "'flux_integral' vector postprocessor.\n\n"
230 : "Length of reporter: ",
231 1 : moose_flux.size(),
232 : "\n",
233 : "Length of 'boundary': ",
234 1 : _boundary->size());
235 :
236 549 : for (std::size_t b = 0; b < _boundary->size(); ++b)
237 : {
238 592 : _console << "[boundary " << Moose::stringify((*_boundary)[b])
239 296 : << "]: Normalizing NekRS flux of "
240 592 : << Moose::stringify(nek_flux_sidesets[b] * nek_flux_print_mult)
241 592 : << " to the conserved MOOSE value of " << Moose::stringify(moose_flux[b])
242 296 : << std::endl;
243 :
244 296 : checkInitialFluxValues(nek_flux_sidesets[b], moose_flux[b]);
245 : }
246 :
247 253 : total_moose_flux = std::accumulate(moose_flux.begin(), moose_flux.end(), 0.0);
248 :
249 : // For the sake of printing diagnostics to the screen regarding the flux normalization,
250 : // we first scale the nek flux by any unit changes and then by the reference flux.
251 : successful_normalization =
252 253 : normalizeFluxBySideset(moose_flux, nek_flux_sidesets, normalized_nek_flux);
253 : }
254 : else
255 : {
256 13060 : auto moose_flux = *_flux_integral;
257 : const double nek_flux =
258 13060 : std::accumulate(nek_flux_sidesets.begin(), nek_flux_sidesets.end(), 0.0);
259 :
260 26120 : _console << "[boundary " << Moose::stringify(*_boundary)
261 13060 : << "]: Normalizing total NekRS flux of "
262 26120 : << Moose::stringify(nek_flux * nek_flux_print_mult)
263 13060 : << " to the conserved MOOSE value of " << Moose::stringify(moose_flux) << std::endl;
264 13060 : _console << Moose::stringify(nek_flux) << " " << nek_flux_print_mult << std::endl;
265 :
266 13060 : checkInitialFluxValues(nek_flux, moose_flux);
267 :
268 13060 : total_moose_flux = moose_flux;
269 :
270 : // For the sake of printing diagnostics to the screen regarding the flux normalization,
271 13060 : successful_normalization = normalizeFlux(moose_flux, nek_flux, normalized_nek_flux);
272 : }
273 :
274 13313 : if (!successful_normalization)
275 1 : mooseError(
276 : "Flux normalization process failed! NekRS integrated flux: ",
277 : normalized_nek_flux,
278 : " MOOSE integrated flux: ",
279 : total_moose_flux,
280 1 : normalizationHint(),
281 : "\n\n- If you set 'conserve_flux_by_sideset = true' and nodes are SHARED by boundaries "
282 : "(like on corners between sidesets), you will end up renormalizing those shared nodes "
283 : "once per sideset that they lie on. There is no guarantee that the total imposed flux "
284 : "would be preserved.");
285 13312 : }
286 :
287 : void
288 13356 : NekBoundaryFlux::checkInitialFluxValues(const Real & nek_flux, const Real & moose_flux) const
289 : {
290 13356 : const Real scale_squared = _nek_mesh->scaling() * _nek_mesh->scaling();
291 13356 : const double nek_flux_print_mult = scale_squared * nekrs::nondimensionalDivisor(field::flux);
292 :
293 : // If before normalization, there is a large difference between the nekRS imposed flux
294 : // and the MOOSE flux, this could mean that there is a poor match between the domains,
295 : // even if neither value is zero. For instance, if you forgot that the nekRS mesh is in
296 : // units of centimeters, but you're coupling to an app based in meters, the fluxes will
297 : // be very different from one another.
298 13356 : if (moose_flux && (std::abs(nek_flux * nek_flux_print_mult - moose_flux) / moose_flux) > 0.25)
299 72 : mooseDoOnce(mooseWarning(
300 : "NekRS flux differs from MOOSE flux by more than 25\%! This is NOT necessarily a problem - "
301 : "but it could indicate that your geometries don't line up properly or something is amiss "
302 : "with your transfer. We recommend opening the output files to visually inspect the flux in "
303 : "both the main and sub applications to check that the fields look correct."));
304 13356 : }
305 :
306 : bool
307 253 : NekBoundaryFlux::normalizeFluxBySideset(const std::vector<double> & moose_integral,
308 : std::vector<double> & nek_integral,
309 : double & normalized_nek_integral)
310 : {
311 : // scale the nek flux to dimensional form for the sake of normalizing against
312 : // a dimensional MOOSE flux
313 549 : for (auto & i : nek_integral)
314 296 : i *= _reference_flux_integral;
315 :
316 253 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
317 253 : mesh_t * mesh = nekrs::temperatureMesh();
318 253 : auto nek_boundary_coupling = _nek_mesh->boundaryCoupling();
319 :
320 146405 : for (int k = 0; k < nek_boundary_coupling.total_n_faces; ++k)
321 : {
322 146152 : if (nek_boundary_coupling.process[k] == nekrs::commRank())
323 : {
324 25264 : int i = nek_boundary_coupling.element[k];
325 25264 : int j = nek_boundary_coupling.face[k];
326 :
327 25264 : int face_id = mesh->EToB[i * mesh->Nfaces + j];
328 25264 : auto it = std::find(_boundary->begin(), _boundary->end(), face_id);
329 : auto b_index = it - _boundary->begin();
330 :
331 : // avoid divide-by-zero
332 : double ratio = 1.0;
333 25264 : if (std::abs(nek_integral[b_index]) > _abs_tol)
334 24512 : ratio = moose_integral[b_index] / nek_integral[b_index];
335 :
336 25264 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
337 :
338 670768 : for (int v = 0; v < mesh->Nfp; ++v)
339 : {
340 645504 : int id = mesh->vmapM[offset + v];
341 645504 : nrs->usrwrk[indices.flux + id] *= ratio;
342 : }
343 : }
344 : }
345 :
346 : // check that the normalization worked properly - confirm against dimensional form
347 253 : auto integrals = nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
348 253 : normalized_nek_integral =
349 253 : std::accumulate(integrals.begin(), integrals.end(), 0.0) * _reference_flux_integral;
350 : double total_moose_integral = std::accumulate(moose_integral.begin(), moose_integral.end(), 0.0);
351 : bool low_rel_err =
352 253 : std::abs(total_moose_integral) > _abs_tol
353 253 : ? std::abs(normalized_nek_integral - total_moose_integral) / total_moose_integral <
354 249 : _rel_tol
355 : : true;
356 253 : bool low_abs_err = std::abs(normalized_nek_integral - total_moose_integral) < _abs_tol;
357 :
358 506 : return low_rel_err || low_abs_err;
359 253 : }
360 :
361 : bool
362 13060 : NekBoundaryFlux::normalizeFlux(const double moose_integral,
363 : double nek_integral,
364 : double & normalized_nek_integral)
365 : {
366 : // scale the nek flux to dimensional form for the sake of normalizing against
367 : // a dimensional MOOSE flux
368 13060 : nek_integral *= _reference_flux_integral;
369 :
370 13060 : auto nek_boundary_coupling = _nek_mesh->boundaryCoupling();
371 :
372 : // avoid divide-by-zero
373 13060 : if (std::abs(nek_integral) < _abs_tol)
374 : return true;
375 :
376 12876 : nrs_t * nrs = (nrs_t *)nekrs::nrsPtr();
377 12876 : mesh_t * mesh = nekrs::temperatureMesh();
378 :
379 12876 : const double ratio = moose_integral / nek_integral;
380 :
381 2049276 : for (int k = 0; k < nek_boundary_coupling.total_n_faces; ++k)
382 : {
383 2036400 : if (nek_boundary_coupling.process[k] == nekrs::commRank())
384 : {
385 390624 : int i = nek_boundary_coupling.element[k];
386 390624 : int j = nek_boundary_coupling.face[k];
387 390624 : int offset = i * mesh->Nfaces * mesh->Nfp + j * mesh->Nfp;
388 :
389 7235064 : for (int v = 0; v < mesh->Nfp; ++v)
390 : {
391 6844440 : int id = mesh->vmapM[offset + v];
392 6844440 : nrs->usrwrk[indices.flux + id] *= ratio;
393 : }
394 : }
395 : }
396 :
397 : // check that the normalization worked properly - confirm against dimensional form
398 12876 : auto integrals = nekrs::usrwrkSideIntegral(indices.flux, *_boundary, nek_mesh::all);
399 12876 : normalized_nek_integral =
400 12876 : std::accumulate(integrals.begin(), integrals.end(), 0.0) * _reference_flux_integral;
401 12876 : bool low_rel_err = std::abs(normalized_nek_integral - moose_integral) / moose_integral < _rel_tol;
402 12876 : bool low_abs_err = std::abs(normalized_nek_integral - moose_integral) < _abs_tol;
403 :
404 12876 : return low_rel_err || low_abs_err;
405 13060 : }
406 :
407 : #endif
|
