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-2025 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
28namespace oomph
29{
30
31////////////////////////////////////////////////////////////////////////
32////////////////////////////////////////////////////////////////////////
33////////////////////////////////////////////////////////////////////////
34
35
36
37//=====================================================================
38/// Sarah's boundary layer solution for flow in oscillating ring
39//=====================================================================
40namespace SarahBL
41{
42
43 double epsilon,alpha,A,Omega,N;
44
45/* The options were : operatorarrow */
46#include <math.h>
47double 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>
90double 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*
122Omega*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>
130double 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>
153double 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*
175t16)/8.0);
176 }
177}
178
179/* The options were : operatorarrow */
180#include <math.h>
181double 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>
214double 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+
274t12*N*t19*t22*t14*t25)*t35)-epsilon*t12*(t14*(Omega*t21*t47*t49-t52)+t14*(t55*
275t21*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*
277t25))*epsilon);
278 }
279}
280
281/* The options were : operatorarrow */
282#include <math.h>
283double 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-
343t12*N*t19*t22*t14*t25)*t35)+epsilon*t12*(t14*(Omega*t21*t47*t49-t52)+t14*(t55*
344t21*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>
352double 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>
373double 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
394 /// Residual function for buckled ring
395void 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
417 /// Exact solution: x,y,u,v,p
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
479 /// Full exact solution: x,y,u,v,p,du/dt,dv/dt,diss
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
oomph::SarahBL::X_sarah
double X_sarah(double rho, double zeta, double t)
Definition osc_ring_sarah_asymptotics.h:352
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::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::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