Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AdjointAdvection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: AdjointAdvection.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: Evaluation of the adjoint advective term
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
42 "Adjoint", AdjointAdvection::create);
43
44/**
45 *
46 */
48{
49}
50
52 const int nConvectiveFields,
54 const Array<OneD, Array<OneD, NekDouble>> &advVel,
55 const Array<OneD, Array<OneD, NekDouble>> &inarray,
56 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time,
57 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
58 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
59{
60 ASSERTL1(nConvectiveFields == inarray.size(),
61 "Number of convective fields and Inarray are not compatible");
62
63 int nPointsTot = fields[0]->GetNpoints();
64 int ndim = advVel.size();
65 int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
66 int nDerivs = (m_halfMode) ? 2 : m_spacedim;
67
69 int nScalar = nConvectiveFields - ndim;
71
72 for (int i = 0; i < ndim; ++i)
73 {
74 if (fields[i]->GetWaveSpace() && !m_singleMode && !m_halfMode)
75 {
76 velocity[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
77 fields[i]->HomogeneousBwdTrans(nPointsTot, advVel[i], velocity[i]);
78 }
79 else
80 {
81 velocity[i] = advVel[i];
82 }
83 }
84 if (nScalar > 0) // add for temperature field
85 {
86 for (int jj = ndim; jj < nConvectiveFields; ++jj)
87 {
88 scalar[jj - ndim] = inarray[jj];
89 }
90 }
91
93 for (int i = 0; i < nDerivs; ++i)
94 {
95 grad[i] = Array<OneD, NekDouble>(nPointsTot);
96 }
97
98 // Evaluation of the base flow for periodic cases
99 if (m_slices > 1)
100 {
101 for (int i = 0; i < ndim; ++i)
102 {
103 UpdateBase(m_interp[i], m_baseflow[i], m_period - time);
104 UpdateGradBase(i, fields[i]);
105 }
106 }
107
108 // Evaluate the linearised advection term
109 for (int i = 0; i < nConvectiveFields; ++i)
110 {
111 // Calculate gradient
112 switch (nDerivs)
113 {
114 case 1:
115 {
116 fields[i]->PhysDeriv(inarray[i], grad[0]);
117 }
118 break;
119 case 2:
120 {
121 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
122 }
123 break;
124 case 3:
125 {
126 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1], grad[2]);
127 if (m_multipleModes)
128 {
129 // transform gradients into physical Fourier space
130 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[0],
131 grad[0]);
132 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[1],
133 grad[1]);
134 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[2],
135 grad[2]);
136 }
137 }
138 break;
139 }
140
141 // Momentum field advection
142 if (i < ndim)
143 {
144 // Calculate -U_j du'_i/dx_j
145 Vmath::Vmul(nPointsTot, grad[0], 1, m_baseflow[0], 1, outarray[i],
146 1);
147 for (int j = 1; j < nDerivs; ++j)
148 {
149 Vmath::Vvtvp(nPointsTot, grad[j], 1, m_baseflow[j], 1,
150 outarray[i], 1, outarray[i], 1);
151 }
152 Vmath::Neg(nPointsTot, outarray[i], 1);
153
154 // Add u'_j U_j/ dx_i
155 int lim = (m_halfMode) ? 2 : ndim;
156 if ((m_halfMode || m_singleMode) && i == 2)
157 {
158 lim = 0;
159 }
160 for (int j = 0; j < lim; ++j)
161 {
162 Vmath::Vvtvp(nPointsTot, m_gradBase[j * nBaseDerivs + i], 1,
163 velocity[j], 1, outarray[i], 1, outarray[i], 1);
164 }
165 // Add Tprime*Grad_Tbase in u, v equations
166 if (nScalar > 0 && i < ndim)
167 {
168 for (int s = 0; s < nScalar; ++s)
169 {
170 Vmath::Vvtvp(nPointsTot,
171 m_gradBase[(ndim + s) * nBaseDerivs + i], 1,
172 scalar[s], 1, outarray[i], 1, outarray[i], 1);
173 }
174 }
175 }
176 // Scalar Field Advection
177 else
178 {
179 // Calculate -U_j du'_i/dx_j
180 Vmath::Vmul(nPointsTot, grad[0], 1, m_baseflow[0], 1, outarray[i],
181 1);
182 for (int j = 1; j < nDerivs; ++j)
183 {
184 Vmath::Vvtvp(nPointsTot, grad[j], 1, m_baseflow[j], 1,
185 outarray[i], 1, outarray[i], 1);
186 }
187 Vmath::Neg(nPointsTot, outarray[i], 1);
188 }
189
190 if (m_multipleModes)
191 {
192 fields[i]->HomogeneousFwdTrans(nPointsTot, outarray[i],
193 outarray[i]);
194 }
195 Vmath::Neg(nPointsTot, outarray[i], 1);
196 }
197}
198
199} // namespace Nektar
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Advects a vector field.
static std::string className
Name of class.
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void UpdateBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble time)
NekDouble m_period
period length
void UpdateGradBase(const int var, const MultiRegions::ExpListSharedPtr &field)
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
Array< OneD, Array< OneD, NekDouble > > m_gradBase
bool m_multipleModes
flag to determine if use multiple mode or not
bool m_halfMode
flag to determine if use half mode or not
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366