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