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 "NekMeshDeformation.h"
22 : #include "DisplacedProblem.h"
23 :
24 : registerMooseObject("CardinalApp", NekMeshDeformation);
25 :
26 : extern nekrs::usrwrkIndices indices;
27 :
28 : InputParameters
29 86 : NekMeshDeformation::validParams()
30 : {
31 86 : auto params = FieldTransferBase::validParams();
32 86 : params.addClassDescription("Reads/writes mesh deformation between NekRS and MOOSE.");
33 86 : return params;
34 0 : }
35 :
36 45 : NekMeshDeformation::NekMeshDeformation(const InputParameters & parameters)
37 45 : : FieldTransferBase(parameters)
38 : {
39 45 : if (_direction == "to_nek" /* && nekrs::hasBlendingSolver() */)
40 : {
41 45 : if (_usrwrk_slot.size() != 3)
42 0 : paramError("usrwrk_slot", "For mesh deformation, 'usrwrk_slot' must be of length 3");
43 :
44 45 : addExternalVariable(_usrwrk_slot[0], _variable + "_x");
45 45 : addExternalVariable(_usrwrk_slot[1], _variable + "_y");
46 90 : addExternalVariable(_usrwrk_slot[2], _variable + "_z");
47 :
48 45 : indices.mesh_velocity_x = _usrwrk_slot[0] * nekrs::fieldOffset();
49 45 : indices.mesh_velocity_y = _usrwrk_slot[1] * nekrs::fieldOffset();
50 45 : indices.mesh_velocity_z = _usrwrk_slot[2] * nekrs::fieldOffset();
51 : }
52 : else
53 : {
54 : // technically, we do not need usrwrk_slot when not using the blending mesh solver,
55 : // since we directly apply the displacements to the entire mesh volume. However,
56 : // the base class is set up to require usrwrk_slot whenever we have displacements
57 : // being applied. Because the mesh deformation is not a widely used feature yet,
58 : // we ignore this annoyance. If someone wants to fix in future with a redesign
59 : // they are welcome to.
60 : // checkUnusedParam(parameters, "usrwrk_slot", "not using the blending mesh solver");
61 : }
62 :
63 45 : if (_direction == "from_nek")
64 0 : paramError("direction",
65 : "The NekMeshDeformation currently only supports transfers 'to_nek'; contact the "
66 : "Cardinal developer team if you require reading of NekRS mesh coordinates.");
67 :
68 45 : if (_app.actionWarehouse().displacedMesh() && !nekrs::hasMovingMesh())
69 1 : mooseWarning(
70 1 : "Your NekRSMesh has 'displacements', but '" + _nek_problem.casename() +
71 : ".par' does not have a\n"
72 : "solver in the [MESH] block! The displacements transferred to NekRS will be unused.");
73 :
74 44 : if (nekrs::hasMovingMesh() && _nek_mesh->exactMirror())
75 0 : mooseError("An exact mesh mirror is not yet implemented for the boundary mesh solver.");
76 :
77 44 : int n_per_surf = _nek_mesh->exactMirror() ? std::pow(_nek_mesh->nekNumQuadraturePoints1D(), 2.0)
78 44 : : _nek_mesh->numVerticesPerSurface();
79 44 : int n_per_vol = _nek_mesh->exactMirror() ? std::pow(_nek_mesh->nekNumQuadraturePoints1D(), 3.0)
80 44 : : _nek_mesh->numVerticesPerVolume();
81 :
82 44 : if (nekrs::hasMovingMesh())
83 : {
84 44 : int n_entries = _nek_mesh->volume() ? n_per_vol : n_per_surf;
85 44 : _displacement_x = (double *)calloc(n_entries, sizeof(double));
86 44 : _displacement_y = (double *)calloc(n_entries, sizeof(double));
87 44 : _displacement_z = (double *)calloc(n_entries, sizeof(double));
88 :
89 44 : if (nekrs::hasBlendingSolver())
90 18 : _mesh_velocity_elem = (double *)calloc(n_entries, sizeof(double));
91 : }
92 :
93 44 : auto boundary = _nek_mesh->boundary();
94 44 : if (!boundary && nekrs::hasBlendingSolver())
95 2 : mooseError("'" + _nek_problem.casename() +
96 : ".par' has a solver in the [MESH] block. This solver uses\n"
97 : "boundary displacement values from MOOSE to move the NekRS mesh. Please indicate\n"
98 : "the 'boundary' for which mesh motion is coupled from MOOSE to NekRS.");
99 :
100 43 : if (nekrs::hasBlendingSolver())
101 : {
102 : bool has_one_mv_bc = false;
103 18 : for (const auto & b : *boundary)
104 : {
105 17 : if (nekrs::isMovingMeshBoundary(b))
106 : {
107 : has_one_mv_bc = true;
108 : break;
109 : }
110 : }
111 :
112 17 : if (!has_one_mv_bc)
113 1 : mooseError("For boundary-coupled moving mesh problems, you need at least one "
114 1 : "boundary in '" +
115 1 : _nek_problem.casename() +
116 : ".par'\nto be of the type 'codedFixedValue'"
117 : " in the [MESH] block.");
118 : }
119 :
120 42 : if (!_nek_mesh->volume() && nekrs::hasUserMeshSolver())
121 1 : mooseError(
122 1 : "'" + _nek_problem.casename() +
123 : ".par' has 'solver = user' in the [MESH] block. With this solver,\n"
124 : "displacement values are sent to every GLL point in NekRS's volume. If you only are "
125 : "building\n"
126 : "a boundary mesh mirror, it's possible that some displacement values could result\n"
127 : "in negative Jacobians if a sideset moves beyond the bounds of an undeformed element.\n"
128 : "To eliminate this possibility, please enable 'volume = true' for NekRSMesh and send a\n"
129 : "whole-domain displacement to NekRS.");
130 :
131 41 : if (nekrs::hasBlendingSolver())
132 16 : _nek_mesh->initializePreviousDisplacements();
133 :
134 41 : if (nekrs::hasUserMeshSolver())
135 25 : _nek_mesh->saveInitialVolMesh();
136 41 : }
137 :
138 80 : NekMeshDeformation::~NekMeshDeformation()
139 : {
140 40 : freePointer(_displacement_x);
141 40 : freePointer(_displacement_y);
142 40 : freePointer(_displacement_z);
143 40 : freePointer(_mesh_velocity_elem);
144 80 : }
145 :
146 : void
147 904 : NekMeshDeformation::sendDataToNek()
148 : {
149 904 : if (nekrs::hasUserMeshSolver())
150 104 : sendVolumeDeformationToNek();
151 800 : else if (nekrs::hasBlendingSolver())
152 800 : sendBoundaryDeformationToNek();
153 : else
154 0 : mooseError("Unhandled mesh solver case in NekMeshDeformation!");
155 :
156 904 : _nek_problem.getDisplacedProblem()->updateMesh();
157 904 : }
158 :
159 : void
160 104 : NekMeshDeformation::sendVolumeDeformationToNek()
161 : {
162 104 : _console << "Sending volume deformation to NekRS" << std::endl;
163 :
164 104 : auto d = nekrs::nondimensionalDivisor(field::x_displacement);
165 104 : auto a = nekrs::nondimensionalAdditive(field::x_displacement);
166 6760 : for (unsigned int e = 0; e < _nek_mesh->numVolumeElems(); e++)
167 : {
168 : // We can only write into the nekRS scratch space if that face is "owned" by the current process
169 6656 : if (nekrs::commRank() != _nek_mesh->volumeCoupling().processor_id(e))
170 5632 : continue;
171 :
172 1024 : _nek_problem.mapVolumeDataToNekVolume(
173 1024 : e, _variable_number[_variable + "_x"], d, a, &_displacement_x);
174 1024 : _nek_problem.writeVolumeSolution(
175 1024 : e, field::x_displacement, _displacement_x, &(_nek_mesh->nek_initial_x()));
176 :
177 1024 : _nek_problem.mapVolumeDataToNekVolume(
178 1024 : e, _variable_number[_variable + "_y"], d, a, &_displacement_y);
179 1024 : _nek_problem.writeVolumeSolution(
180 1024 : e, field::y_displacement, _displacement_y, &(_nek_mesh->nek_initial_y()));
181 :
182 1024 : _nek_problem.mapVolumeDataToNekVolume(
183 1024 : e, _variable_number[_variable + "_z"], d, a, &_displacement_z);
184 1024 : _nek_problem.writeVolumeSolution(
185 1024 : e, field::z_displacement, _displacement_z, &(_nek_mesh->nek_initial_z()));
186 : }
187 104 : }
188 :
189 : void
190 800 : NekMeshDeformation::sendBoundaryDeformationToNek()
191 : {
192 800 : _console << "Sending boundary deformation to NekRS..." << std::endl;
193 :
194 800 : auto d = nekrs::nondimensionalDivisor(field::x_displacement);
195 800 : auto a = nekrs::nondimensionalAdditive(field::x_displacement);
196 800 : if (!_nek_mesh->volume())
197 : {
198 179600 : for (unsigned int e = 0; e < _nek_mesh->numSurfaceElems(); e++)
199 : {
200 : // We can only write into the nekRS scratch space if that face is "owned" by the current
201 : // process
202 179200 : if (nekrs::commRank() != _nek_mesh->boundaryCoupling().processor_id(e))
203 134400 : continue;
204 :
205 44800 : _nek_problem.mapFaceDataToNekFace(
206 44800 : e, _variable_number[_variable + "_x"], d, a, &_displacement_x);
207 44800 : calculateMeshVelocity(e, field::mesh_velocity_x);
208 44800 : _nek_problem.writeBoundarySolution(e, field::mesh_velocity_x, _mesh_velocity_elem);
209 :
210 44800 : _nek_problem.mapFaceDataToNekFace(
211 44800 : e, _variable_number[_variable + "_y"], d, a, &_displacement_y);
212 44800 : calculateMeshVelocity(e, field::mesh_velocity_y);
213 44800 : _nek_problem.writeBoundarySolution(e, field::mesh_velocity_y, _mesh_velocity_elem);
214 :
215 44800 : _nek_problem.mapFaceDataToNekFace(
216 44800 : e, _variable_number[_variable + "_z"], d, a, &_displacement_z);
217 44800 : calculateMeshVelocity(e, field::mesh_velocity_z);
218 44800 : _nek_problem.writeBoundarySolution(e, field::mesh_velocity_z, _mesh_velocity_elem);
219 : }
220 : }
221 : else
222 : {
223 538000 : for (unsigned int e = 0; e < _nek_mesh->numVolumeElems(); ++e)
224 : {
225 : // We can only write into the nekRS scratch space if that face is "owned" by the current
226 : // process
227 537600 : if (nekrs::commRank() != _nek_mesh->volumeCoupling().processor_id(e))
228 403200 : continue;
229 :
230 134400 : _nek_problem.mapFaceDataToNekVolume(
231 134400 : e, _variable_number[_variable + "_x"], d, a, &_displacement_x);
232 134400 : calculateMeshVelocity(e, field::mesh_velocity_x);
233 134400 : _nek_problem.writeVolumeSolution(e, field::mesh_velocity_x, _mesh_velocity_elem);
234 :
235 134400 : _nek_problem.mapFaceDataToNekVolume(
236 134400 : e, _variable_number[_variable + "_y"], d, a, &_displacement_y);
237 134400 : calculateMeshVelocity(e, field::mesh_velocity_y);
238 134400 : _nek_problem.writeVolumeSolution(e, field::mesh_velocity_y, _mesh_velocity_elem);
239 :
240 134400 : _nek_problem.mapFaceDataToNekVolume(
241 134400 : e, _variable_number[_variable + "_z"], d, a, &_displacement_z);
242 134400 : calculateMeshVelocity(e, field::mesh_velocity_z);
243 134400 : _nek_problem.writeVolumeSolution(e, field::mesh_velocity_z, _mesh_velocity_elem);
244 : }
245 : }
246 800 : }
247 :
248 : void
249 537600 : NekMeshDeformation::calculateMeshVelocity(int e, const field::NekWriteEnum & field)
250 : {
251 : int len =
252 537600 : _nek_mesh->volume() ? _nek_mesh->numVerticesPerVolume() : _nek_mesh->numVerticesPerSurface();
253 :
254 537600 : double dt = _nek_problem.transientExecutioner()->getTimeStepper()->getCurrentDT();
255 :
256 : double *displacement = nullptr, *prev_disp = nullptr;
257 : field::NekWriteEnum disp_field;
258 :
259 537600 : switch (field)
260 : {
261 179200 : case field::mesh_velocity_x:
262 179200 : displacement = _displacement_x;
263 : prev_disp = _nek_mesh->prev_disp_x().data();
264 : disp_field = field::x_displacement;
265 179200 : break;
266 179200 : case field::mesh_velocity_y:
267 179200 : displacement = _displacement_y;
268 : prev_disp = _nek_mesh->prev_disp_y().data();
269 : disp_field = field::y_displacement;
270 179200 : break;
271 179200 : case field::mesh_velocity_z:
272 179200 : displacement = _displacement_z;
273 : prev_disp = _nek_mesh->prev_disp_z().data();
274 : disp_field = field::z_displacement;
275 179200 : break;
276 0 : default:
277 0 : mooseError("Unhandled NekWriteEnum in NekRSProblem::calculateMeshVelocity!\n");
278 : }
279 :
280 537600 : auto reference_v = nekrs::nondimensionalDivisor(field::velocity);
281 12633600 : for (int i = 0; i < len; i++)
282 12096000 : _mesh_velocity_elem[i] = (displacement[i] - prev_disp[(e * len) + i]) / dt / reference_v;
283 :
284 537600 : _nek_mesh->updateDisplacement(e, displacement, disp_field);
285 537600 : }
286 :
287 : #endif
|
