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