osc_ring_sarah_asymptotics.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 
27 
28 namespace oomph
29 {
30 
34 
35 
36 
37 //=====================================================================
39 //=====================================================================
40 namespace SarahBL
41 {
42 
43  double epsilon,alpha,A,Omega,N;
44 
45 /* The options were : operatorarrow */
46 #include <math.h>
47 double Diss_sarah(double rho,double zeta,double t)
48 {
49  double t1;
50  double t10;
51  double t11;
52  double t13;
53  double t14;
54  double t19;
55  double t2;
56  double t20;
57  double t21;
58  double t22;
59  double t25;
60  double t26;
61  double t39;
62  double t43;
63  double t5;
64  double t6;
65  double t8;
66  {
67  t1 = alpha*alpha;
68  t2 = epsilon*epsilon;
69  t5 = sin(N*zeta);
70  t6 = t5*t5;
71  t8 = Omega*Omega;
72  t10 = 1.0+A;
73  t11 = t10*t10;
74  t13 = sqrt(2.0);
75  t14 = sqrt(Omega);
76  t19 = t13*t14*alpha*(1.0-rho)/2.0;
77  t20 = exp(-t19);
78  t21 = t20*t20;
79  t22 = Omega*t;
80  t25 = cos(-t19+t22+0.3141592653589793E1/4.0);
81  t26 = t25*t25;
82  t39 = cos(-t19+t22);
83  t43 = cos(t22);
84  return(t1*t2*t6*t8*Omega*t11*t21*t26+4.0*alpha*t2*t6*t14*Omega*t10*t20*t25*
85 (N-1.0)*(Omega*t10*t20*t39/2.0-Omega*t43));
86  }
87 }
88 /* The options were : operatorarrow */
89 #include <math.h>
90 double Kin_energy_sarah(double t)
91 {
92  double t1;
93  double t11;
94  double t13;
95  double t16;
96  double t18;
97  double t2;
98  double t20;
99  double t27;
100  double t28;
101  double t31;
102  double t33;
103  double t5;
104  double t6;
105  double t7;
106  {
107  t1 = epsilon*epsilon;
108  t2 = Omega*Omega;
109  t5 = Omega*t;
110  t6 = cos(t5);
111  t7 = t6*t6;
112  t11 = 1/alpha;
113  t13 = sqrt(Omega);
114  t16 = 1.0+A;
115  t18 = 0.3141592653589793E1/4.0;
116  t20 = sin(t5+t18);
117  t27 = t16*t16;
118  t28 = 1/t13;
119  t31 = sin(2.0*t5+t18);
120  t33 = sqrt(2.0);
121  return(t1*(0.3141592653589793E1*t2/N*t7/2.0+t11*0.3141592653589793E1*t13*
122 Omega*t16*t6*t20)+0.3141592653589793E1*t1*t2*t11*(t27*(t28*t31/4.0+t33*t28/4.0
123 )-2.0*t28*t16*t6*t20)/2.0);
124  }
125 }
126 
127 
128 /* The options were : operatorarrow */
129 #include <math.h>
130 double P_sarah(double rho,double zeta,double t)
131 {
132  double t1;
133  double t12;
134  double t19;
135  double t4;
136  double t6;
137  double t8;
138  double t9;
139  {
140  t1 = Omega*Omega;
141  t4 = pow(rho,1.0*N);
142  t6 = cos(N*zeta);
143  t8 = Omega*t;
144  t9 = sin(t8);
145  t12 = sqrt(Omega);
146  t19 = cos(t8+0.3141592653589793E1/4.0);
147  return(epsilon*(t1/N*t4*t6*t9-t12*Omega*(1.0+A)*t4*t6*t19/alpha));
148  }
149 }
150 
151 /* The options were : operatorarrow */
152 #include <math.h>
153 double Total_Diss_lead_sarah(double t)
154 {
155  double t1;
156  double t10;
157  double t11;
158  double t12;
159  double t13;
160  double t16;
161  double t4;
162  double t5;
163  double t6;
164  {
165  t1 = epsilon*epsilon;
166  t4 = sqrt(2.0);
167  t5 = Omega*Omega;
168  t6 = sqrt(Omega);
169  t10 = pow(1.0+A,2.0);
170  t11 = Omega*t;
171  t12 = sin(t11);
172  t13 = cos(t11);
173  t16 = t13*t13;
174  return(-alpha*t1*0.3141592653589793E1*t4*t6*t5*t10*(-1.0+2.0*t12*t13-2.0*
175 t16)/8.0);
176  }
177 }
178 
179 /* The options were : operatorarrow */
180 #include <math.h>
181 double Total_Diss_sarah(double t)
182 {
183  double t1;
184  double t14;
185  double t15;
186  double t18;
187  double t19;
188  double t2;
189  double t20;
190  double t4;
191  double t6;
192  double t7;
193  double t8;
194  {
195  t1 = Omega*Omega;
196  t2 = 0.3141592653589793E1*t1;
197  t4 = epsilon*epsilon;
198  t6 = Omega*t;
199  t7 = cos(t6);
200  t8 = t7*t7;
201  t14 = sqrt(2.0);
202  t15 = sqrt(Omega);
203  t18 = 1.0+A;
204  t19 = t18*t18;
205  t20 = sin(t6);
206  return(4.0*t2*(N-1.0)*t4*t8-alpha*t4*0.3141592653589793E1*t14*t15*t1*t19*(
207 -1.0+2.0*t20*t7-2.0*t8)/8.0+t4*t18*t2*(-10.0*A*t8-A+8.0*A*N*t8+22.0*t8-1.0-24.0
208 *N*t8)/8.0);
209  }
210 }
211 
212 /* The options were : operatorarrow */
213 #include <math.h>
214 double U_sarah(double rho,double zeta,double t)
215 {
216  double t1;
217  double t12;
218  double t14;
219  double t19;
220  double t21;
221  double t22;
222  double t23;
223  double t25;
224  double t3;
225  double t35;
226  double t4;
227  double t41;
228  double t43;
229  double t46;
230  double t47;
231  double t49;
232  double t5;
233  double t52;
234  double t55;
235  double t57;
236  double t59;
237  double t6;
238  double t62;
239  double t64;
240  double t8;
241  double t87;
242  double t89;
243  double t9;
244  {
245  t1 = cos(zeta);
246  t3 = N-1.0;
247  t4 = pow(rho,1.0*t3);
248  t5 = N*zeta;
249  t6 = cos(t5);
250  t8 = Omega*t;
251  t9 = cos(t8);
252  t12 = sin(zeta);
253  t14 = sin(t5);
254  t19 = sqrt(Omega);
255  t21 = 1.0+A;
256  t22 = t21*t4;
257  t23 = 0.3141592653589793E1/4.0;
258  t25 = sin(t8+t23);
259  t35 = 1/alpha;
260  t41 = sqrt(2.0);
261  t43 = 1.0-rho;
262  t46 = t41*t19*alpha*t43/2.0;
263  t47 = exp(-t46);
264  t49 = cos(-t46+t8);
265  t52 = Omega*t9;
266  t55 = N*t19;
267  t57 = sin(-t46+t8+t23);
268  t59 = t47*t57-t25;
269  t62 = Omega*alpha;
270  t64 = -t43*t3;
271  t87 = t9*(1.0+t64);
272  t89 = N*t35;
273  return(epsilon*(t1*Omega*t4*t6*t9+t12*Omega*t4*t14*t9+(t1*N*t19*t22*t6*t25+
274 t12*N*t19*t22*t14*t25)*t35)-epsilon*t12*(t14*(Omega*t21*t47*t49-t52)+t14*(t55*
275 t21*t59-t62*t64*t9)*t35)+epsilon*t1*(t52*t6+t6*(-t55*t21*t59-t62*t43*t3*t9)*t35
276 )-(-t12*(-Omega*t14*t87-t89*t19*t21*t14*t25)+t1*(Omega*t6*t87+t89*t6*t19*t21*
277 t25))*epsilon);
278  }
279 }
280 
281 /* The options were : operatorarrow */
282 #include <math.h>
283 double V_sarah(double rho,double zeta,double t)
284 {
285  double t1;
286  double t12;
287  double t14;
288  double t19;
289  double t21;
290  double t22;
291  double t23;
292  double t25;
293  double t3;
294  double t35;
295  double t4;
296  double t41;
297  double t43;
298  double t46;
299  double t47;
300  double t49;
301  double t5;
302  double t52;
303  double t55;
304  double t57;
305  double t59;
306  double t6;
307  double t62;
308  double t64;
309  double t8;
310  double t87;
311  double t89;
312  double t9;
313  {
314  t1 = sin(zeta);
315  t3 = N-1.0;
316  t4 = pow(rho,1.0*t3);
317  t5 = N*zeta;
318  t6 = cos(t5);
319  t8 = Omega*t;
320  t9 = cos(t8);
321  t12 = cos(zeta);
322  t14 = sin(t5);
323  t19 = sqrt(Omega);
324  t21 = 1.0+A;
325  t22 = t21*t4;
326  t23 = 0.3141592653589793E1/4.0;
327  t25 = sin(t8+t23);
328  t35 = 1/alpha;
329  t41 = sqrt(2.0);
330  t43 = 1.0-rho;
331  t46 = t41*t19*alpha*t43/2.0;
332  t47 = exp(-t46);
333  t49 = cos(-t46+t8);
334  t52 = Omega*t9;
335  t55 = N*t19;
336  t57 = sin(-t46+t8+t23);
337  t59 = t47*t57-t25;
338  t62 = Omega*alpha;
339  t64 = -t43*t3;
340  t87 = t9*(1.0+t64);
341  t89 = N*t35;
342  return(epsilon*(t1*Omega*t4*t6*t9-t12*Omega*t4*t14*t9+(t1*N*t19*t22*t6*t25-
343 t12*N*t19*t22*t14*t25)*t35)+epsilon*t12*(t14*(Omega*t21*t47*t49-t52)+t14*(t55*
344 t21*t59-t62*t64*t9)*t35)+epsilon*t1*(t52*t6+t6*(-t55*t21*t59-t62*t43*t3*t9)*t35
345 )-(t12*(-Omega*t14*t87-t89*t19*t21*t14*t25)+t1*(Omega*t6*t87+t89*t6*t19*t21*t25
346 ))*epsilon);
347  }
348 }
349 
350 /* The options were : operatorarrow */
351 #include <math.h>
352 double X_sarah(double rho,double zeta,double t)
353 {
354  double t1;
355  double t10;
356  double t3;
357  double t5;
358  double t6;
359  double t8;
360  {
361  t1 = cos(zeta);
362  t3 = sin(Omega*t);
363  t5 = N*zeta;
364  t6 = cos(t5);
365  t8 = sin(t5);
366  t10 = sin(zeta);
367  return(rho*(t1+epsilon*t3*(t6*t1-A*t8*t10)));
368  }
369 }
370 
371 /* The options were : operatorarrow */
372 #include <math.h>
373 double Y_sarah(double rho,double zeta,double t)
374 {
375  double t1;
376  double t10;
377  double t3;
378  double t5;
379  double t6;
380  double t8;
381  {
382  t1 = sin(zeta);
383  t3 = sin(Omega*t);
384  t5 = N*zeta;
385  t6 = cos(t5);
386  t8 = sin(t5);
387  t10 = cos(zeta);
388  return(rho*(t1+epsilon*t3*(t6*t1+A*t8*t10)));
389  }
390 }
391 
392 
393 
395 void buckled_ring_residual(const Vector<double>& params,
396  const Vector<double>& unknowns,
397  Vector<double>& residuals)
398 {
399 
400  // Parameters
401  double t=params[0];
402  double xx=params[1];
403  double yy=params[2];
404 
405  // Unknowns
406  double rho=unknowns[0];
407  double zeta=unknowns[1];
408 
409  // Residuals
410  residuals[0]=X_sarah(rho,zeta,t)-xx;
411  residuals[1]=Y_sarah(rho,zeta,t)-yy;
412 
413  }
414 
415 
416 
418  void exact_soln(const double& time,
419  const Vector<double>& x,
420  Vector<double>& soln)
421  {
422 
423  // Primary variables only
424  soln.resize(3);
425 
426  // Get rho and zeta with black-box Newton method:
427 
428  // Parameters: time, x, y
429  Vector<double> params(3);
430  params[0]=time;
431  params[1]=x[0];
432  params[2]=x[1];
433 
434 
435  // Unknowns rho, zeta
436  Vector<double> unknowns(2);
437 
438  // Initial guess:
439  double rho=sqrt(x[0]*x[0]+x[1]*x[1]);
440  double zeta=0.5*MathematicalConstants::Pi;
441  if (std::abs(x[0])>1e-4) zeta=atan(x[1]/x[0]);
442 
443  // Copy across
444  unknowns[0]=rho;
445  unknowns[1]=zeta;
446 
447  // Call Newton solver
448  BlackBoxFDNewtonSolver::black_box_fd_newton_solve(
449  buckled_ring_residual,params,unknowns);
450 
451  // Extract rho, zeta
452 
453  // Copy across
454  rho=unknowns[0];
455  zeta=unknowns[1];
456 
457 
458  // Double check the position
459  double dx=X_sarah(rho,zeta,time)-x[0];
460  double dy=Y_sarah(rho,zeta,time)-x[1];
461  double error=sqrt(dx*dx+dy*dy);
462  if (error>1.0e-6)
463  {
464  std::cout << "Trafo error " << error << std::endl;
465  //assert(false);
466  }
467 
468  // Velocities
469  soln[0]=U_sarah(rho,zeta,time);
470  soln[1]=V_sarah(rho,zeta,time);
471 
472  // Pressure on viscous scale
473  soln[2]=P_sarah(rho,zeta,time)*(alpha*alpha);
474 
475  }
476 
477 
478 
480  void full_exact_soln(const double& time,
481  const Vector<double>& x,
482  Vector<double>& soln)
483  {
484 
485 
486  // Full solution
487  soln.resize(6);
488 
489  // Get rho and zeta with black-box Newton method:
490 
491  // Parameters: time, x, y
492  Vector<double> params(3);
493  params[0]=time;
494  params[1]=x[0];
495  params[2]=x[1];
496 
497 
498  // Unknowns rho, zeta
499  Vector<double> unknowns(2);
500 
501  // Initial guess:
502  double rho=sqrt(x[0]*x[0]+x[1]*x[1]);
503  double zeta=0.5*MathematicalConstants::Pi;
504  if (std::abs(x[0])>1e-4) zeta=atan(x[1]/x[0]);
505 
506  // Copy across
507  unknowns[0]=rho;
508  unknowns[1]=zeta;
509 
510  // Call Newton solver
511  BlackBoxFDNewtonSolver::black_box_fd_newton_solve(
512  buckled_ring_residual,params,unknowns);
513 
514  // Extract rho, zeta
515 
516  // Copy across
517  rho=unknowns[0];
518  zeta=unknowns[1];
519 
520 
521  // Double check the position
522  double dx=X_sarah(rho,zeta,time)-x[0];
523  double dy=Y_sarah(rho,zeta,time)-x[1];
524  double error=sqrt(dx*dx+dy*dy);
525  if (error>1.0e-6)
526  {
527  std::cout << "Trafo error " << error << std::endl;
528  //assert(false);
529  }
530 
531  // Velocities
532  soln[0]=U_sarah(rho,zeta,time);
533  soln[1]=V_sarah(rho,zeta,time);
534 
535  // Pressure on viscous scale
536  soln[2]=P_sarah(rho,zeta,time)*(alpha*alpha);
537 
538  // du/dt dv/dt not coded up yet...
539  soln[3]=0.0; //DUDT(xx,y,time);
540  soln[4]=0.0; //DVDT(xx,y,time);
541 
542  if (rho<1.0e-6)
543  {
544  rho=1e-6;
545  }
546  soln[5]=Diss_sarah(rho,zeta,time);
547 
548 
549  }
550 
551 }
552 
553 
554 }
555 
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
Eigen::bfloat16_impl::atan
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 atan(const bfloat16 &a)
Definition: BFloat16.h:636
Eigen::bfloat16_impl::exp
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 exp(const bfloat16 &a)
Definition: BFloat16.h:615
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
Particles2023AnalysisHung.params
dictionary params
Definition: Particles2023AnalysisHung.py:35
SarahBL
Sarah's boundary layer solution for flow in oscillating ring.
Definition: navier_stokes/osc_ring/osc_ring_alg.cc:72
calibrate.error
int error
Definition: calibrate.py:297
oomph::BlackBoxFDNewtonSolver::black_box_fd_newton_solve
void black_box_fd_newton_solve(ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &unknowns, JacobianFctPt jacobian_fct)
Definition: black_box_newton_solver.cc:67
oomph::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition: oomph_utilities.h:157
oomph::SarahBL::X_sarah
double X_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:352
oomph::SarahBL::epsilon
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43
oomph::SarahBL::Omega
double Omega
Definition: osc_ring_sarah_asymptotics.h:43
oomph::SarahBL::Diss_sarah
double Diss_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:47
oomph::SarahBL::Kin_energy_sarah
double Kin_energy_sarah(double t)
Definition: osc_ring_sarah_asymptotics.h:90
oomph::SarahBL::N
double N
Definition: osc_ring_sarah_asymptotics.h:43
oomph::SarahBL::A
double A
Definition: osc_ring_sarah_asymptotics.h:43
oomph::SarahBL::U_sarah
double U_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:214
oomph::SarahBL::full_exact_soln
void full_exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Full exact solution: x,y,u,v,p,du/dt,dv/dt,diss.
Definition: osc_ring_sarah_asymptotics.h:480
oomph::SarahBL::P_sarah
double P_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:130
oomph::SarahBL::buckled_ring_residual
void buckled_ring_residual(const Vector< double > &params, const Vector< double > &unknowns, Vector< double > &residuals)
Residual function for buckled ring.
Definition: osc_ring_sarah_asymptotics.h:395
oomph::SarahBL::Total_Diss_lead_sarah
double Total_Diss_lead_sarah(double t)
Definition: osc_ring_sarah_asymptotics.h:153
oomph::SarahBL::Y_sarah
double Y_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:373
oomph::SarahBL::exact_soln
void exact_soln(const double &time, const Vector< double > &x, Vector< double > &soln)
Exact solution: x,y,u,v,p.
Definition: osc_ring_sarah_asymptotics.h:418
oomph::SarahBL::alpha
double alpha
Definition: osc_ring_sarah_asymptotics.h:43
oomph::SarahBL::Total_Diss_sarah
double Total_Diss_sarah(double t)
Definition: osc_ring_sarah_asymptotics.h:181
oomph::SarahBL::V_sarah
double V_sarah(double rho, double zeta, double t)
Definition: osc_ring_sarah_asymptotics.h:283
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition: AnisotropicHookean.h:10
plotDoE.x
list x
Definition: plotDoE.py:28
plotPSD.t
t
Definition: plotPSD.py:36