Nektar++
ALEHelper.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ALEHelper.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Helper class for ALE process
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include "ALEHelper.h"
38
39namespace Nektar::SolverUtils
40{
41
43 [[maybe_unused]] int spaceDim,
44 [[maybe_unused]] Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
45{
46}
47
49 [[maybe_unused]] int spaceDim,
50 [[maybe_unused]] Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
51{
52 // Create ALE objects for each interface zone
53 if (fields[0]->GetGraph() !=
54 nullptr) // homogeneous graphs are missing the graph data
55 {
56 for (auto &zone : fields[0]->GetGraph()->GetMovement()->GetZones())
57 {
58 switch (zone.second->GetMovementType())
59 {
61 m_ALEs.emplace_back(ALEFixedShPtr(
63 zone.second)));
64 break;
66 m_ALEs.emplace_back(ALETranslateShPtr(
68 zone.second)));
69 m_ALESolver = true;
70 break;
72 m_ALEs.emplace_back(ALERotateShPtr(
74 zone.second)));
75 m_ALESolver = true;
76 break;
78 WARNINGL0(false,
79 "Zone cannot have movement type of 'None'.");
80 default:
81 break;
82 }
83 }
84 }
85
86 // Update grid velocity
88}
89
91{
92
93 // Reset grid velocity to 0
94 for (int i = 0; i < m_spaceDim; ++i)
95 {
96 std::fill(m_gridVelocity[i].begin(), m_gridVelocity[i].end(), 0.0);
97 }
98 if (m_ALESolver)
99 {
100 // Now update for each movement zone, adding the grid velocities
101 for (auto &ALE : m_ALEs)
102 {
103 ALE->UpdateGridVel(time, m_fieldsALE, m_gridVelocity);
104 }
105 }
106}
107
110{
111 if (m_ALESolver)
112 {
113 const int nm = m_fieldsALE[0]->GetNcoeffs();
115
116 // Premultiply each field by the mass matrix
117 for (int i = 0; i < m_fieldsALE.size(); ++i)
118 {
119 fields[i] = Array<OneD, NekDouble>(nm);
120 m_fieldsALE[i]->GeneralMatrixOp(mkey, m_fieldsALE[i]->GetCoeffs(),
121 fields[i]);
122 }
123 }
124}
125
126/**
127* @brief Update m_fields with u^n by multiplying by inverse mass
128 matrix. That's then used in e.g. checkpoint output and L^2 error
129 calculation.
130*/
132 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &traceNormals,
134 [[maybe_unused]] NekDouble time)
135{
136 // @TODO: Look at geometric factor and junk only what is needed
137 // @TODO: Look at collections and see if they offer a speed up
138 for (int i = 0; i < m_fieldsALE.size(); ++i)
139 {
140 m_fieldsALE[i]->MultiplyByElmtInvMass(
141 fields[i],
142 m_fieldsALE[i]->UpdateCoeffs()); // @TODO: Potentially matrix free?
143 m_fieldsALE[i]->BwdTrans(m_fieldsALE[i]->GetCoeffs(),
144 m_fieldsALE[i]->UpdatePhys());
145 }
146}
147
149 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
151{
152 const int nc = m_fieldsALE[0]->GetNcoeffs();
153 int nVariables = inarray.size();
154
155 // General idea is that we are time-integrating the quantity (Mu), so we
156 // need to multiply input by inverse mass matrix to get coefficients u,
157 // and then backwards transform to physical space so we can apply the DG
158 // operator.
160 for (int i = 0; i < nVariables; ++i)
161 {
162 outarray[i] = Array<OneD, NekDouble>(m_fieldsALE[0]->GetNpoints());
163 m_fieldsALE[i]->MultiplyByElmtInvMass(inarray[i], tmp);
164 m_fieldsALE[i]->BwdTrans(tmp, outarray[i]);
165 }
166}
167
169 Array<OneD, Array<OneD, NekDouble>> &traceNormals)
170{
171 // Only move if timestepped
172 if (time == m_prevStageTime)
173 {
174 return;
175 }
176
177 auto curvedEdges = m_fieldsALE[0]->GetGraph()->GetCurvedEdges();
178 auto curvedFaces = m_fieldsALE[0]->GetGraph()->GetCurvedFaces();
179
181 timer.Start();
182 m_fieldsALE[0]->GetGraph()->GetMovement()->PerformMovement(
183 time); // @TODO: Moved out of loop!
184 timer.Stop();
185 timer.AccumulateRegion("Movement::PerformMovement");
186
187 // The order of the resets below is v important to avoid errors
188 for (auto &field : m_fieldsALE)
189 {
190 field->ResetMatrices();
191 }
192
193 // Loop over all elements and faces and edges and reset geometry
194 // information. Only need to do this on the first field as the geometry
195 // information is shared.
196 for (auto &zone : m_fieldsALE[0]->GetGraph()->GetMovement()->GetZones())
197 {
198 if (zone.second->GetMoved())
199 {
200 auto conEl = zone.second->GetConstituentElements();
201 for (const auto &i : conEl)
202 {
203 for (const auto &j : i)
204 {
205 j->ResetNonRecursive(curvedEdges, curvedFaces);
206 }
207 }
208
209 // We need to rebuild geometric factors on the trace elements
210 for (const auto &i : conEl[m_fieldsALE[0]->GetShapeDimension() -
211 1]) // This only takes the trace elements
212 {
213 m_fieldsALE[0]
214 ->GetTrace()
215 ->GetExpFromGeomId(i->GetGlobalID())
216 ->Reset();
217 }
218 }
219 }
220
221 for (auto &field : m_fieldsALE)
222 {
223 for (auto &zone : field->GetGraph()->GetMovement()->GetZones())
224 {
225 if (zone.second->GetMoved())
226 {
227 auto conEl = zone.second->GetConstituentElements();
228 // Loop over zone elements expansions and rebuild geometric
229 // factors
230 for (const auto &i :
231 conEl[0]) // This only takes highest dimensioned elements
232 {
233 field->GetExpFromGeomId(i->GetGlobalID())->Reset();
234 }
235 }
236 }
237 }
238
239 for (auto &zone : m_fieldsALE[0]->GetGraph()->GetMovement()->GetZones())
240 {
241 if (zone.second->GetMoved())
242 {
243 auto conEl = zone.second->GetConstituentElements();
244 // Loop over zone elements expansions and rebuild geometric factors
245 // and recalc trace normals
246 for (const auto &i :
247 conEl[0]) // This only takes highest dimensioned elements
248 {
249 int nfaces = m_fieldsALE[0]
250 ->GetExpFromGeomId(i->GetGlobalID())
251 ->GetNtraces();
252 for (int j = 0; j < nfaces; ++j)
253 {
254 m_fieldsALE[0]
255 ->GetExpFromGeomId(i->GetGlobalID())
256 ->ComputeTraceNormal(j);
257 }
258 }
259 }
260 }
261
262 for (auto &field : m_fieldsALE)
263 {
264 // Reset collections (despite the default being eNoCollection it does
265 // remember the last auto-tuned values), eNoImpType gives lots of output
266 field->CreateCollections(Collections::eNoCollection);
267 }
268
269 // Reload new trace normals in to the solver cache
270 m_fieldsALE[0]->GetTrace()->GetNormals(traceNormals);
271
272 // Recompute grid velocity.
274
275 // Updates trace grid velocity
276 for (int i = 0; i < m_gridVelocityTrace.size(); ++i)
277 {
278 m_fieldsALE[0]->ExtractTracePhys(m_gridVelocity[i],
280 }
281
282 // Set the flag to exchange coords in InterfaceMapDG to true
283 m_fieldsALE[0]->GetGraph()->GetMovement()->GetCoordExchangeFlag() = true;
284
285 m_prevStageTime = time;
286}
287
290{
291 return m_gridVelocityTrace;
292}
293
295 : m_zone(std::static_pointer_cast<SpatialDomains::ZoneFixed>(zone))
296{
297}
298
300 [[maybe_unused]] NekDouble time,
301 [[maybe_unused]] Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
302 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &gridVelocity)
303{
304}
305
307 : m_zone(std::static_pointer_cast<SpatialDomains::ZoneTranslate>(zone))
308{
309}
310
312 [[maybe_unused]] NekDouble time,
314 Array<OneD, Array<OneD, NekDouble>> &gridVelocity)
315{
316 auto vel = m_zone->GetVel(time);
317 auto exp = fields[0]->GetExp();
318
319 auto elements = m_zone->GetElements();
320 for (auto &el : elements)
321 {
322 int indx = fields[0]->GetElmtToExpId(el->GetGlobalID());
323 int offset = fields[0]->GetPhys_Offset(indx);
324 auto expansion = (*exp)[indx];
325
326 int nq = expansion->GetTotPoints();
327 for (int i = 0; i < nq; ++i)
328 {
329 for (int j = 0; j < gridVelocity.size(); ++j)
330 {
331 gridVelocity[j][offset + i] += vel[j];
332 }
333 }
334 }
335}
336
338 : m_zone(std::static_pointer_cast<SpatialDomains::ZoneRotate>(zone))
339{
340}
341
343 [[maybe_unused]] NekDouble time,
344 [[maybe_unused]] Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
345 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &gridVelocity)
346{
347 auto angVel = m_zone->GetAngularVel(time);
348 auto axis = m_zone->GetAxis();
349 auto origin = m_zone->GetOrigin();
350
351 auto exp = fields[0]->GetExp();
352
353 auto elements = m_zone->GetElements();
354 for (auto &el : elements)
355 {
356 int indx = fields[0]->GetElmtToExpId(el->GetGlobalID());
357 int offset = fields[0]->GetPhys_Offset(indx);
358 auto expansion = (*exp)[indx];
359
360 int nq = expansion->GetTotPoints();
361
362 Array<OneD, NekDouble> xc(nq, 0.0), yc(nq, 0.0), zc(nq, 0.0);
363 Array<OneD, NekDouble> norm(3, 0.0);
364 expansion->GetCoords(xc, yc, zc);
365 for (int i = 0; i < nq; ++i)
366 {
367 // Vector from origin to point
368 NekDouble xpointMinOrigin = xc[i] - origin(0);
369 NekDouble ypointMinOrigin = yc[i] - origin(1);
370 NekDouble zpointMinOrigin = zc[i] - origin(2);
371
372 // Vector orthogonal to plane formed by axis and point
373 // We negate angVel here as by convention a positive angular
374 // velocity is counter-clockwise
375 norm[0] = (ypointMinOrigin * axis[2] - zpointMinOrigin * axis[1]) *
376 (-angVel);
377 norm[1] = (zpointMinOrigin * axis[0] - xpointMinOrigin * axis[2]) *
378 (-angVel);
379 norm[2] = (xpointMinOrigin * axis[1] - ypointMinOrigin * axis[0]) *
380 (-angVel);
381 for (int j = 0; j < gridVelocity.size(); ++j)
382 {
383 gridVelocity[j][offset + i] = norm[j];
384 }
385 }
386 }
387}
388
390 std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
391 std::vector<std::string> &variables)
392
393{
394 int nCoeffs = m_fieldsALE[0]->GetNcoeffs();
395 // Adds extra output variables for grid velocity
396 std::string gridVarName[3] = {"gridVx", "gridVy", "gridVz"};
397 for (int i = 0; i < m_spaceDim; ++i)
398 {
399 Array<OneD, NekDouble> gridVel(nCoeffs, 0.0);
400 m_fieldsALE[0]->FwdTransLocalElmt(m_gridVelocity[i], gridVel);
401 fieldcoeffs.emplace_back(gridVel);
402 variables.emplace_back(gridVarName[i]);
403 }
404}
405
406} // namespace Nektar::SolverUtils
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:70
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Describes a matrix with ordering defined by a local to global map.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fieldsALE
Definition: ALEHelper.h:88
virtual SOLVER_UTILS_EXPORT void v_ALEPreMultiplyMass(Array< OneD, Array< OneD, NekDouble > > &fields)
Definition: ALEHelper.cpp:108
virtual SOLVER_UTILS_EXPORT void v_ALEInitObject(int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ALEHelper.cpp:42
SOLVER_UTILS_EXPORT void ALEDoElmtInvMass(Array< OneD, Array< OneD, NekDouble > > &traceNormals, Array< OneD, Array< OneD, NekDouble > > &fields, NekDouble time)
Update m_fields with u^n by multiplying by inverse mass matrix. That's then used in e....
Definition: ALEHelper.cpp:131
SOLVER_UTILS_EXPORT void InitObject(int spaceDim, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: ALEHelper.cpp:48
Array< OneD, Array< OneD, NekDouble > > m_gridVelocityTrace
Definition: ALEHelper.h:90
SOLVER_UTILS_EXPORT void ALEDoElmtInvMassBwdTrans(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ALEHelper.cpp:148
SOLVER_UTILS_EXPORT void ExtraFldOutputGridVelocity(std::vector< Array< OneD, NekDouble > > &fieldcoeffs, std::vector< std::string > &variables)
Definition: ALEHelper.cpp:389
SOLVER_UTILS_EXPORT void MoveMesh(const NekDouble &time, Array< OneD, Array< OneD, NekDouble > > &traceNormals)
Definition: ALEHelper.cpp:168
virtual SOLVER_UTILS_EXPORT void v_UpdateGridVelocity(const NekDouble &time)
Definition: ALEHelper.cpp:90
SOLVER_UTILS_EXPORT const Array< OneD, const Array< OneD, NekDouble > > & GetGridVelocityTrace()
Definition: ALEHelper.cpp:289
Array< OneD, Array< OneD, NekDouble > > m_gridVelocity
Definition: ALEHelper.h:89
std::vector< ALEBaseShPtr > m_ALEs
Definition: ALEHelper.h:91
const NekPoint< NekDouble > origin
std::shared_ptr< ALETranslate > ALETranslateShPtr
Definition: ALEHelper.h:157
std::shared_ptr< ALERotate > ALERotateShPtr
Definition: ALEHelper.h:158
std::shared_ptr< ALEFixed > ALEFixedShPtr
Definition: ALEHelper.h:156
std::shared_ptr< ZoneBase > ZoneBaseShPtr
Definition: Zones.h:171
double NekDouble
STL namespace.
void v_UpdateGridVel(const NekDouble time, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, Array< OneD, Array< OneD, NekDouble > > &gridVelocity) final
Definition: ALEHelper.cpp:299
ALEFixed(SpatialDomains::ZoneBaseShPtr zone)
Definition: ALEHelper.cpp:294
SpatialDomains::ZoneRotateShPtr m_zone
Definition: ALEHelper.h:153
ALERotate(SpatialDomains::ZoneBaseShPtr zone)
Definition: ALEHelper.cpp:337
void v_UpdateGridVel(const NekDouble time, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, Array< OneD, Array< OneD, NekDouble > > &gridVelocity) final
Definition: ALEHelper.cpp:342
SpatialDomains::ZoneTranslateShPtr m_zone
Definition: ALEHelper.h:140
ALETranslate(SpatialDomains::ZoneBaseShPtr zone)
Definition: ALEHelper.cpp:306
void v_UpdateGridVel(const NekDouble time, Array< OneD, MultiRegions::ExpListSharedPtr > &fields, Array< OneD, Array< OneD, NekDouble > > &gridVelocity) final
Definition: ALEHelper.cpp:311