Functions | Variables
PerturbationSolution Namespace Reference

Christian's perturbation solution. More...

Functions

double get_radial_distance (const Vector< double > &x)
 function to compute radial distance of this point More...
 
double get_azimuthal_angle (const Vector< double > &x)
 function to return angle of radial line from positive x-axis More...
 
void polar_to_cartesian_velocity (const Vector< double > u_polar, const double phi, Vector< double > &u_cart)
 function to convert velocity vectors from polar to Cartesian coords More...
 
void leading_order_veloc_and_pressure (const double r, const double phi, Vector< double > &u_0, double &p_0)
 Leading order solution (u and p are in polar coordinates) More...
 
void first_order_veloc_and_pressure (const double r, const double phi, Vector< double > &u_1, double &p_1)
 First order correction (u and p are in polar coordinates) More...
 
void second_order_veloc_and_pressure (const double r, const double phi, Vector< double > &u_2, double &p_2)
 Second order correction (u and p are in polar coordinates) More...
 
void full_soln (const Vector< double > &x, Vector< double > &u, double &p, const unsigned &nterms=2)
 
void perturbation_soln_for_plot (const Vector< double > &x, Vector< double > &soln)
 Combined solution (u,v,p) for plotting. More...
 
void first_order_perturbation_for_plot (const Vector< double > &x, Vector< double > &soln)
 First order-perturbation for plotting. More...
 
void second_order_perturbation_for_plot (const Vector< double > &x, Vector< double > &soln)
 Second order-perturbation for plotting. More...
 

Variables

void(* veloc_and_pressure_term_pt [3])(const double, const double, Vector< double > &, double &)
 array of function pointers for each term in expansion More...
 
unsigned N_terms_for_plot =2
 Number of terms in perturbation expansion for plot. More...
 

Detailed Description

Christian's perturbation solution.

Function Documentation

◆ first_order_perturbation_for_plot()

void PerturbationSolution::first_order_perturbation_for_plot ( const Vector< double > &  x,
Vector< double > &  soln 
)

First order-perturbation for plotting.

210  {
211  soln.resize(3);
212 
213  // get radial distance to this point from origin
214  double r = get_radial_distance(x);
215 
216  // get azimuthal angle of this point
217  double phi = get_azimuthal_angle(x);
218 
219  // vector to hold each velocity term (Cartesian)
220  Vector<double> u_i(2);
221 
222  // vector to hold velocity term (polar coordinates)
223  Vector<double> u_polar(2);
224 
225  // variable to hold each pressure term
226  double p_i = 0;
227 
228  // compute the term (in polar coordinates)
229  veloc_and_pressure_term_pt[0](r, phi, u_polar, p_i);
230 
231  // convert coordinates to Cartesians
232  polar_to_cartesian_velocity(u_polar, phi, u_i);
233 
234  // multiply by appropriate power of Re and add it to the solution
235  soln[0]=-u_i[0];
236  soln[1]=-u_i[1];
237  soln[2]=-p_i;
238 
239  for (unsigned i=0;i<3;i++)
240  {
241  if (std::isinf(soln[i])) soln[i]=200.0;
242  if (std::isnan(soln[i])) soln[i]=200.0;
243  }
244 
245  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
isnan
#define isnan(X)
Definition: main.h:109
isinf
#define isinf(X)
Definition: main.h:110
PerturbationSolution::get_azimuthal_angle
double get_azimuthal_angle(const Vector< double > &x)
function to return angle of radial line from positive x-axis
Definition: finite_re_perturbation.h:43
PerturbationSolution::polar_to_cartesian_velocity
void polar_to_cartesian_velocity(const Vector< double > u_polar, const double phi, Vector< double > &u_cart)
function to convert velocity vectors from polar to Cartesian coords
Definition: finite_re_perturbation.h:49
PerturbationSolution::veloc_and_pressure_term_pt
void(* veloc_and_pressure_term_pt[3])(const double, const double, Vector< double > &, double &)
array of function pointers for each term in expansion
Definition: finite_re_perturbation.h:136
PerturbationSolution::get_radial_distance
double get_radial_distance(const Vector< double > &x)
function to compute radial distance of this point
Definition: finite_re_perturbation.h:37
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
plotDoE.x
list x
Definition: plotDoE.py:28

References get_azimuthal_angle(), get_radial_distance(), i, isinf, isnan, polar_to_cartesian_velocity(), UniformPSDSelfTest::r, veloc_and_pressure_term_pt, and plotDoE::x.

◆ first_order_veloc_and_pressure()

void PerturbationSolution::first_order_veloc_and_pressure ( const double  r,
const double  phi,
Vector< double > &  u_1,
double &  p_1 
)

First order correction (u and p are in polar coordinates)

75  {
76  // first order velocity term, r component
77  u_1[0] = (r*(-2*pow(Pi,2)*(8 + pow(Pi,2)) +
78  2*(8*pow(phi,2)*(-4 + pow(Pi,2)) - 8*phi*Pi*(2 + pow(Pi,2)) + pow(Pi,2)*(8 + pow(Pi,2)))*
79  cos(2*phi) + (32*phi - 72*Pi - 64*pow(phi,2)*Pi + 24*phi*pow(Pi,2) + 6*pow(Pi,3) + pow(Pi,5))*
80  sin(2*phi)))/(32.*pow(-4 + pow(Pi,2),2));
81 
82  // phi component
83  u_1[1] = (r*(Pi*(24 - 14*pow(Pi,2) - pow(Pi,4) + 4*phi*Pi*(8 + pow(Pi,2))) +
84  (-64*pow(phi,2)*Pi + 8*phi*(12 + pow(Pi,2)) + Pi*(-24 + 14*pow(Pi,2) + pow(Pi,4)))*cos(2*phi) -
85  2*(24 + 10*pow(Pi,2) + pow(Pi,4) + 8*pow(phi,2)*(-4 + pow(Pi,2)) - 8*phi*Pi*(6 + pow(Pi,2)))*
86  sin(2*phi)))/(32.*pow(-4 + pow(Pi,2),2));
87 
88  // first order pressure term
89  p_1 = -(pow(Pi,2)*(-8 + pow(Pi,2))*log(r))/(4.*pow(-4 + pow(Pi,2),2));
90  }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
Eigen::ArrayBase::pow
EIGEN_DEVICE_FUNC const GlobalUnaryPowReturnType< Derived, ScalarExponent > pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
Definition: GlobalFunctions.h:137
BiharmonicTestFunctions2::Pi
double Pi
Definition: two_d_biharmonic.cc:235
Eigen::bfloat16_impl::log
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 log(const bfloat16 &a)
Definition: BFloat16.h:618

References cos(), Eigen::bfloat16_impl::log(), BiharmonicTestFunctions2::Pi, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, and sin().

◆ full_soln()

void PerturbationSolution::full_soln ( const Vector< double > &  x,
Vector< double > &  u,
double &  p,
const unsigned &  nterms = 2 
)

Full solution. Final optional argument specifies number of terms in expansion (Stokes or Stokes plus leading order correction in Re)

153  {
154  // get radial distance to this point from origin
155  double r = get_radial_distance(x);
156 
157  // get azimuthal angle of this point
158  double phi = get_azimuthal_angle(x);
159 
160  // vector to hold each velocity term (Cartesian)
161  Vector<double> u_i(2);
162 
163  // vector to hold velocity term (polar coordinates)
164  Vector<double> u_polar(2);
165 
166  // variable to hold each pressure term
167  double p_i = 0;
168 
169  // initialise solutions
170  u[0] = 0.0;
171  u[1] = 0.0;
172  p = 0.0;
173 
174  // -------------------------------------------
175  // compute each term and add it to the solution
176  for (unsigned i = 0; i < nterms; i++)
177  {
178  // compute the term (in polar coordinates)
179  veloc_and_pressure_term_pt[i](r, phi, u_polar, p_i);
180 
181  // convert coordinates to Cartesians
182  polar_to_cartesian_velocity(u_polar, phi, u_i);
183 
184  // multiply by appropriate power of Re and add it to the solution
185  u[0] -= pow(Re, i) * u_i[0];
186  u[1] -= pow(Re, i) * u_i[1];
187  p -= pow(Re, i) * p_i;
188  }
189  }
p
float * p
Definition: Tutorial_Map_using.cpp:9
GlobalPhysicalVariables::Re
double Re
Reynolds number.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:61

References get_azimuthal_angle(), get_radial_distance(), i, p, polar_to_cartesian_velocity(), Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, GlobalPhysicalVariables::Re, veloc_and_pressure_term_pt, and plotDoE::x.

Referenced by perturbation_soln_for_plot().

◆ get_azimuthal_angle()

double PerturbationSolution::get_azimuthal_angle ( const Vector< double > &  x)

function to return angle of radial line from positive x-axis

44  {
45  return atan2(x[1], x[0]);
46  }
atan2
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:139

References atan2(), and plotDoE::x.

Referenced by first_order_perturbation_for_plot(), full_soln(), and second_order_perturbation_for_plot().

◆ get_radial_distance()

double PerturbationSolution::get_radial_distance ( const Vector< double > &  x)

function to compute radial distance of this point

38  {
39  return sqrt(x[0]*x[0] + x[1]*x[1]);
40  }
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134

References sqrt(), and plotDoE::x.

Referenced by first_order_perturbation_for_plot(), full_soln(), and second_order_perturbation_for_plot().

◆ leading_order_veloc_and_pressure()

void PerturbationSolution::leading_order_veloc_and_pressure ( const double  r,
const double  phi,
Vector< double > &  u_0,
double &  p_0 
)

Leading order solution (u and p are in polar coordinates)

60  {
61  // leading order velocity term, r component
62  u_0[0] =
63  (((4 + 2*phi*Pi - pow(Pi,2))*cos(phi) + 2*(-2*phi + Pi)*sin(phi)))/(-4 + pow(Pi,2));
64 
65  // phi component
66  u_0[1] = ((-4*phi*cos(phi) + Pi*(-2*phi + Pi)*sin(phi)))/(-4 + pow(Pi,2));
67 
68  // leading order pressure term
69  p_0 = (4*(2*cos(phi) + Pi*sin(phi)))/((-4 + pow(Pi,2))*r);
70  }

References cos(), BiharmonicTestFunctions2::Pi, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, and sin().

◆ perturbation_soln_for_plot()

void PerturbationSolution::perturbation_soln_for_plot ( const Vector< double > &  x,
Vector< double > &  soln 
)

Combined solution (u,v,p) for plotting.

198  {
199  soln.resize(3);
200  double p=0.0;
201  full_soln(x,soln,p,N_terms_for_plot);
202  if (std::isinf(p)) p=200.0;
203  if (std::isnan(p)) p=200.0;
204  soln[2]=p;
205  }
PerturbationSolution::N_terms_for_plot
unsigned N_terms_for_plot
Number of terms in perturbation expansion for plot.
Definition: finite_re_perturbation.h:193
PerturbationSolution::full_soln
void full_soln(const Vector< double > &x, Vector< double > &u, double &p, const unsigned &nterms=2)
Definition: finite_re_perturbation.h:150

References full_soln(), isinf, isnan, N_terms_for_plot, p, and plotDoE::x.

◆ polar_to_cartesian_velocity()

void PerturbationSolution::polar_to_cartesian_velocity ( const Vector< double >  u_polar,
const double  phi,
Vector< double > &  u_cart 
)

function to convert velocity vectors from polar to Cartesian coords

52  {
53  u_cart[0] = u_polar[0] * cos(phi) - u_polar[1] * sin(phi);
54  u_cart[1] = u_polar[0] * sin(phi) + u_polar[1] * cos(phi);
55  }

References cos(), and sin().

Referenced by first_order_perturbation_for_plot(), full_soln(), and second_order_perturbation_for_plot().

◆ second_order_perturbation_for_plot()

void PerturbationSolution::second_order_perturbation_for_plot ( const Vector< double > &  x,
Vector< double > &  soln 
)

Second order-perturbation for plotting.

250  {
251  soln.resize(3);
252 
253  // get radial distance to this point from origin
254  double r = get_radial_distance(x);
255 
256  // get azimuthal angle of this point
257  double phi = get_azimuthal_angle(x);
258 
259  // vector to hold each velocity term (Cartesian)
260  Vector<double> u_i(2);
261 
262  // vector to hold velocity term (polar coordinates)
263  Vector<double> u_polar(2);
264 
265  // variable to hold each pressure term
266  double p_i = 0;
267 
268  // compute the term (in polar coordinates)
269  veloc_and_pressure_term_pt[1](r, phi, u_polar, p_i);
270 
271  // convert coordinates to Cartesians
272  polar_to_cartesian_velocity(u_polar, phi, u_i);
273 
274  // multiply by appropriate power of Re and add it to the solution
275  soln[0]=-u_i[0];
276  soln[1]=-u_i[1];
277  soln[2]=-p_i;
278 
279  for (unsigned i=0;i<3;i++)
280  {
281  if (std::isinf(soln[i])) soln[i]=200.0;
282  if (std::isnan(soln[i])) soln[i]=200.0;
283  }
284 
285  }

References get_azimuthal_angle(), get_radial_distance(), i, isinf, isnan, polar_to_cartesian_velocity(), UniformPSDSelfTest::r, veloc_and_pressure_term_pt, and plotDoE::x.

◆ second_order_veloc_and_pressure()

void PerturbationSolution::second_order_veloc_and_pressure ( const double  r,
const double  phi,
Vector< double > &  u_2,
double &  p_2 
)

Second order correction (u and p are in polar coordinates)

95  {
96  // second order velocity term, r component
97  u_2[0] = (pow(r,2)*(2*(-3952 + 3672*pow(Pi,2) - 51*pow(Pi,4) - 21*pow(Pi,6) +
98  192*pow(phi,3)*Pi*(4 + pow(Pi,2)) + 72*phi*Pi*(2 + pow(Pi,2))*(-8 + 3*pow(Pi,2)) +
99  144*pow(phi,2)*(20 + 5*pow(Pi,2) - 2*pow(Pi,4)))*cos(phi) +
100  2*(3952 - 3672*pow(Pi,2) + 51*pow(Pi,4) + 21*pow(Pi,6) + 288*pow(phi,3)*Pi*(-12 + pow(Pi,2)) -
101  432*pow(phi,2)*(-12 + 5*pow(Pi,2) + pow(Pi,4)) + 24*phi*Pi*(276 + 85*pow(Pi,2) + 9*pow(Pi,4)))*
102  cos(3*phi) + (-768*pow(phi,3)*(4 + pow(Pi,2)) + 144*pow(phi,2)*Pi*(20 + 13*pow(Pi,2)) -
103  36*phi*(-32 + 32*pow(Pi,2) + 42*pow(Pi,4) + pow(Pi,6)) +
104  Pi*(-3648 + 2576*pow(Pi,2) + 564*pow(Pi,4) + 9*pow(Pi,6)))*sin(phi) +
105  (3072*Pi + 1152*pow(phi,3)*(4 - 3*pow(Pi,2)) + 432*pow(phi,2)*Pi*(36 + 5*pow(Pi,2)) -
106  pow(Pi,3)*(2560 + 252*pow(Pi,2) + 27*pow(Pi,4)) +
107  12*phi*(-1024 + 72*pow(Pi,2) + 54*pow(Pi,4) + 3*pow(Pi,6)))*sin(3*phi)))/
108  (4608.*pow(-4 + pow(Pi,2),3));
109 
110  // phi component
111  u_2[1] = (pow(r,2)*((-768*pow(phi,3)*(4 + pow(Pi,2)) + 144*pow(phi,2)*Pi*(-12 + 5*pow(Pi,2)) -
112  36*phi*(-224 - 16*pow(Pi,2) + 10*pow(Pi,4) + pow(Pi,6)) +
113  Pi*(2112 + 1424*pow(Pi,2) + 132*pow(Pi,4) + 9*pow(Pi,6)))*cos(phi) +
114  (-384*pow(phi,3)*(-4 + 3*pow(Pi,2)) + 48*pow(phi,2)*Pi*(156 + 11*pow(Pi,2)) -
115  Pi*(2112 + 1424*pow(Pi,2) + 132*pow(Pi,4) + 9*pow(Pi,6)) +
116  4*phi*(-1856 + 3*pow(Pi,2)*(6 + pow(Pi,2))*(28 + pow(Pi,2))))*cos(3*phi) +
117  2*(-80 - 3960*pow(Pi,2) + 231*pow(Pi,4) + 39*pow(Pi,6) - 192*pow(phi,3)*Pi*(4 + pow(Pi,2)) +
118  144*pow(phi,2)*(12 + 3*pow(Pi,2) + 2*pow(Pi,4)) - 72*phi*Pi*(-40 + 8*pow(Pi,2) + 3*pow(Pi,4)))*
119  sin(phi) - 2*(80 - 888*pow(Pi,2) + 85*pow(Pi,4) + 9*pow(Pi,6) +
120  96*pow(phi,3)*Pi*(-12 + pow(Pi,2)) + 8*phi*Pi*(588 + 107*pow(Pi,2) + 9*pow(Pi,4)) -
121  48*pow(phi,2)*(-52 + 3*pow(Pi,2)*(9 + pow(Pi,2))))*sin(3*phi)))/(1536.*pow(-4 + pow(Pi,2),3));
122 
123  // second order pressure term
124  p_2 = (r*(-((9632 - 6264*pow(Pi,2) + 552*pow(Pi,4) + 33*pow(Pi,6) + 144*phi*Pi*(28 - 13*pow(Pi,2)) +
125  192*pow(phi,3)*Pi*(4 + pow(Pi,2)) - 288*pow(phi,2)*(20 + pow(Pi,2)))*cos(phi)) +
126  9*(-352 + 80*pow(Pi,2) + 22*pow(Pi,4) + pow(Pi,6) + 32*pow(phi,2)*(4 - 3*pow(Pi,2)) +
127  16*phi*Pi*(36 + pow(Pi,2)))*cos(3*phi) +
128  (384*pow(phi,3)*(4 + pow(Pi,2)) - 144*pow(phi,2)*Pi*(-20 + 3*pow(Pi,2) + pow(Pi,4)) +
129  36*phi*(-256 - 72*pow(Pi,2) + 10*pow(Pi,4) + pow(Pi,6)) +
130  Pi*(7152 - 124*pow(Pi,2) + 24*pow(Pi,4) + 9*pow(Pi,6)))*sin(phi) -
131  36*(96*phi - 48*(-2 + pow(phi,2))*Pi - 56*phi*pow(Pi,2) + 4*(4 + pow(phi,2))*pow(Pi,3) -
132  4*phi*pow(Pi,4) + pow(Pi,5))*sin(3*phi)))/(576.*pow(-4 + pow(Pi,2),3));
133  }

References cos(), BiharmonicTestFunctions2::Pi, Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, and sin().

Variable Documentation

◆ N_terms_for_plot

unsigned PerturbationSolution::N_terms_for_plot =2

Number of terms in perturbation expansion for plot.

Referenced by perturbation_soln_for_plot().

◆ veloc_and_pressure_term_pt

void(* PerturbationSolution::veloc_and_pressure_term_pt[3])(const double, const double, Vector< double > &, double &) ( const double  ,
const double  ,
Vector< double > &  ,
double &   
)
Initial value:
=
{ & leading_order_veloc_and_pressure,
& first_order_veloc_and_pressure,
& second_order_veloc_and_pressure
}
PerturbationSolution::second_order_veloc_and_pressure
void second_order_veloc_and_pressure(const double r, const double phi, Vector< double > &u_2, double &p_2)
Second order correction (u and p are in polar coordinates)
Definition: finite_re_perturbation.h:93
PerturbationSolution::leading_order_veloc_and_pressure
void leading_order_veloc_and_pressure(const double r, const double phi, Vector< double > &u_0, double &p_0)
Leading order solution (u and p are in polar coordinates)
Definition: finite_re_perturbation.h:58
PerturbationSolution::first_order_veloc_and_pressure
void first_order_veloc_and_pressure(const double r, const double phi, Vector< double > &u_1, double &p_1)
First order correction (u and p are in polar coordinates)
Definition: finite_re_perturbation.h:73

array of function pointers for each term in expansion

Referenced by first_order_perturbation_for_plot(), full_soln(), and second_order_perturbation_for_plot().