Nektar++
TemplatePressureArea.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TemplatePressureArea.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// License for the specific language governing rights and limitations under
14// Permission is hereby granted, free of charge, to any person obtaining a
15// copy of this software and associated documentation files (the "Software"),
16// to deal in the Software without restriction, including without limitation
17// the rights to use, copy, modify, merge, publish, distribute, sublicense,
18// and/or sell copies of the Software, and to permit persons to whom the
19// Software is furnished to do so, subject to the following conditions:
20//
21// The above copyright notice and this permission notice shall be included
22// in all copies or substantial portions of the Software.
23//
24// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30// DEALINGS IN THE SOFTWARE.
31//
32// Description: TemplatePressureArea class
33//
34///////////////////////////////////////////////////////////////////////////////
35
37
38namespace Nektar
39{
40
44 "Template to add a new Pressure-Area relation");
45
49 : PulseWavePressureArea(pVessel, pSession)
50{
51}
52
54 const NekDouble &A,
55 const NekDouble &A0,
56 const NekDouble &dAUdx,
57 const NekDouble &gamma,
58 const NekDouble &alpha)
59{
60 // Relation here
61
62 /*
63 NOTE: to add a new relation, you also need to add the name of the .h and
64 .cpp files to the CMakeLists.txt file in the PulseWaveSolver folder.
65 */
66}
67
69 const NekDouble A, const NekDouble A0,
70 const NekDouble alpha)
71{
72 // PWV here
73}
74
76 const NekDouble beta, const NekDouble A,
77 const NekDouble A0, const NekDouble alpha)
78{
79 NekDouble I = 0.0;
80 GetCharIntegral(I, beta, A, A0);
81
82 W1 = u + I; // Elastic and assumes u0 = 0
83}
84
86 const NekDouble beta, const NekDouble A,
87 const NekDouble A0, const NekDouble alpha)
88{
89 NekDouble I = 0.0;
90 GetCharIntegral(I, beta, A, A0);
91
92 W2 = u - I; // Elastic and assumes u0 = 0
93}
94
96 const NekDouble W2,
97 const NekDouble beta,
98 const NekDouble A0,
99 const NekDouble alpha)
100{
101 /*
102 If an anayltical formula exists for this, then it should be used here. If
103 not, a numerical method is used.
104 */
105
106 // Iterative Newton-Raphson
107 NekDouble I = 0.0;
108 NekDouble c = 0.0;
109 NekDouble FA = 0.0;
110 NekDouble dFdA = 0.0;
111 NekDouble dA = 0.0;
112
113 int proceed = 1;
114 int iter = 0;
115 int MAX_ITER = 200;
116 NekDouble Tol = 1.0E-10;
117
118 while ((proceed) && (iter < MAX_ITER))
119 {
120 iter += 1;
121
122 GetCharIntegral(I, beta, A, A0, alpha);
123 GetC(c, beta, A, A0, alpha);
124
125 FA = I - ((W1 - W2) / 2);
126 dFdA = c / A;
127
128 dA = FA / dFdA;
129 A -= dA;
130
131 if (sqrt(dA * dA) < Tol)
132 {
133 proceed = 0;
134 }
135 }
136}
137
139 const NekDouble W2)
140{
141 u = (W1 + W2) / 2; // Necessarily the case for all tube laws
142}
143
145 const NekDouble A,
146 const NekDouble A0,
147 const NekDouble alpha)
148{
149 /*
150 If an anayltical formula exists for this, then it should be used here. If
151 not, a numerical method is used.
152 */
153
154 // Numerical integration via the trapezium rule
155 // f(x) = c / A, x1 = A0, x2 = A
156
157 int n = 500;
158 NekDouble h = (A - A0) / n;
159 NekDouble A_n = A0;
160 NekDouble c_n = 0.0;
161
162 for (int i = 1; i < n; ++i)
163 {
164 A_n += h;
165 GetC(c_n, beta, A_n, A0, alpha);
166 I += c_n / A_n;
167 }
168
169 NekDouble c = 0.0;
170 NekDouble c0 = 0.0;
171 GetC(c, beta, A, A0, alpha);
172 GetC(c0, beta, A0, A0, alpha);
173
174 I += ((c / A) + (c0 / A0)) / 2;
175 I *= h;
176}
177
182 const std::string type)
183{
184 /*
185 This is a general method that will work for all tube laws. Simulations will
186 run faster if specific analytical formulae are derived as in the case of the
187 beta law, but it's not necessary.
188 */
189
190 // General formulation
191 if (type == "Bifurcation")
192 {
193 NekMatrix<NekDouble> J(6, 6);
195
196 for (int i = 0; i < 3; ++i)
197 {
198 GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
199 }
200
201 J.SetValue(0, 0, 1);
202 J.SetValue(0, 1, 0);
203 J.SetValue(0, 2, 0);
204 J.SetValue(0, 3, c[0] / Au[0]);
205 J.SetValue(0, 4, 0);
206 J.SetValue(0, 5, 0);
207
208 J.SetValue(1, 0, 0);
209 J.SetValue(1, 1, 1);
210 J.SetValue(1, 2, 0);
211 J.SetValue(1, 3, 0);
212 J.SetValue(1, 4, -c[1] / Au[1]);
213 J.SetValue(1, 5, 0);
214
215 J.SetValue(2, 0, 0);
216 J.SetValue(2, 1, 0);
217 J.SetValue(2, 2, 1);
218 J.SetValue(2, 3, 0);
219 J.SetValue(2, 4, 0);
220 J.SetValue(2, 5, -c[2] / Au[2]);
221
222 J.SetValue(3, 0, Au[0]);
223 J.SetValue(3, 1, -Au[1]);
224 J.SetValue(3, 2, -Au[2]);
225 J.SetValue(3, 3, uu[0]);
226 J.SetValue(3, 4, -uu[1]);
227 J.SetValue(3, 5, -uu[2]);
228
229 J.SetValue(4, 0, 2 * uu[0]);
230 J.SetValue(4, 1, -2 * uu[1]);
231 J.SetValue(4, 2, 0);
232 J.SetValue(4, 3, 2 * c[0] * c[0] / Au[0]);
233 J.SetValue(4, 4, -2 * c[1] * c[1] / Au[1]);
234 J.SetValue(4, 5, 0);
235
236 J.SetValue(5, 0, 2 * uu[0]);
237 J.SetValue(5, 1, 0);
238 J.SetValue(5, 2, -2 * uu[2]);
239 J.SetValue(5, 3, 2 * c[0] * c[0] / Au[0]);
240 J.SetValue(5, 4, 0);
241 J.SetValue(5, 5, -2 * c[2] * c[2] / Au[2]);
242
243 invJ = J;
244 invJ.Invert();
245 }
246 else if (type == "Merge")
247 {
248 NekMatrix<NekDouble> J(6, 6);
250
251 for (int i = 0; i < 3; ++i)
252 {
253 GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
254 }
255
256 J.SetValue(0, 0, 1);
257 J.SetValue(0, 1, 0);
258 J.SetValue(0, 2, 0);
259 J.SetValue(0, 3, -c[0] / Au[0]);
260 J.SetValue(0, 4, 0);
261 J.SetValue(0, 5, 0);
262
263 J.SetValue(1, 0, 0);
264 J.SetValue(1, 1, 1);
265 J.SetValue(1, 2, 0);
266 J.SetValue(1, 3, 0);
267 J.SetValue(1, 4, c[1] / Au[1]);
268 J.SetValue(1, 5, 0);
269
270 J.SetValue(2, 0, 0);
271 J.SetValue(2, 1, 0);
272 J.SetValue(2, 2, 1);
273 J.SetValue(2, 3, 0);
274 J.SetValue(2, 4, 0);
275 J.SetValue(2, 5, c[2] / Au[2]);
276
277 J.SetValue(3, 0, Au[0]);
278 J.SetValue(3, 1, -Au[1]);
279 J.SetValue(3, 2, -Au[2]);
280 J.SetValue(3, 3, uu[0]);
281 J.SetValue(3, 4, -uu[1]);
282 J.SetValue(3, 5, -uu[2]);
283
284 J.SetValue(4, 0, 2 * uu[0]);
285 J.SetValue(4, 1, -2 * uu[1]);
286 J.SetValue(4, 2, 0);
287 J.SetValue(4, 3, 2 * c[0] * c[0] / Au[0]);
288 J.SetValue(4, 4, -2 * c[1] * c[1] / Au[1]);
289 J.SetValue(4, 5, 0);
290
291 J.SetValue(5, 0, 2 * uu[0]);
292 J.SetValue(5, 1, 0);
293 J.SetValue(5, 2, -2 * uu[2]);
294 J.SetValue(5, 3, 2 * c[0] * c[0] / Au[0]);
295 J.SetValue(5, 4, 0);
296 J.SetValue(5, 5, -2 * c[2] * c[2] / Au[2]);
297
298 invJ = J;
299 invJ.Invert();
300 }
301 else if (type == "Interface")
302 {
303 NekMatrix<NekDouble> J(4, 4);
305
306 for (int i = 0; i < 2; ++i)
307 {
308 GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
309 }
310
311 J.SetValue(0, 0, 1);
312 J.SetValue(0, 1, 0);
313 J.SetValue(0, 2, c[0] / Au[0]);
314 J.SetValue(0, 3, 0);
315
316 J.SetValue(1, 0, 0);
317 J.SetValue(1, 1, 1);
318 J.SetValue(1, 2, 0);
319 J.SetValue(1, 3, -c[1] / Au[1]);
320
321 J.SetValue(2, 0, Au[0]);
322 J.SetValue(2, 1, -Au[1]);
323 J.SetValue(2, 2, uu[0]);
324 J.SetValue(2, 3, -uu[1]);
325
326 J.SetValue(3, 0, 2 * uu[0]);
327 J.SetValue(3, 1, -2 * uu[1]);
328 J.SetValue(3, 2, 2 * c[0] * c[0] / Au[0]);
329 J.SetValue(3, 3, -2 * c[1] * c[1] / Au[1]);
330
331 invJ = J;
332 invJ.Invert();
333 }
334}
335
336} // namespace Nektar
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void GetCharIntegral(NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void GetC(NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void v_GetW1(NekDouble &W1, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
virtual void v_GetUFromChars(NekDouble &u, const NekDouble &W1, const NekDouble &W2) override
virtual void v_GetC(NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
virtual void v_GetPressure(NekDouble &P, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &dAUdx, const NekDouble &gamma=0, const NekDouble &alpha=0.5) override
virtual void v_GetAFromChars(NekDouble &A, const NekDouble &W1, const NekDouble &W2, const NekDouble &beta, const NekDouble &A0, const NekDouble &alpha=0.5) override
virtual void v_GetW2(NekDouble &W2, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
virtual void v_GetCharIntegral(NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5) override
static PulseWavePressureAreaSharedPtr create(Array< OneD, MultiRegions::ExpListSharedPtr > &pVessel, const LibUtilities::SessionReaderSharedPtr &pSession)
virtual void v_GetJacobianInverse(NekMatrix< NekDouble > &invJ, const Array< OneD, NekDouble > &Au, const Array< OneD, NekDouble > &uu, const Array< OneD, NekDouble > &beta, const Array< OneD, NekDouble > &A0, const Array< OneD, NekDouble > &alpha, const std::string &type) override
TemplatePressureArea(Array< OneD, MultiRegions::ExpListSharedPtr > pVessel, const LibUtilities::SessionReaderSharedPtr pSession)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
@ P
Monomial polynomials .
Definition: BasisType.h:62
PressureAreaFactory & GetPressureAreaFactory()
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285