cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Autoformatted host_forcing.cc

Browse Source
This commit is contained in:
jpekkila
2019-07-03 17:06:57 +03:00
parent af3a1e211e
commit 945751e585
1 changed files with 55 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files

113
src/standalone/model/host_forcing.cc
View File

@@ -28,17 +28,14 @@
#include "core/math_utils.h" #include "core/math_utils.h"
// The is a wrapper for genering random numbers with a chosen system.
//The is a wrapper for genering random numbers with a chosen system.
AcReal AcReal
get_random_number_01() get_random_number_01()
{ {
//TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ // TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/
return AcReal(rand())/AcReal(RAND_MAX); return AcReal(rand()) / AcReal(RAND_MAX);
} }
AcReal3 AcReal3
cross(const AcReal3& a, const AcReal3& b) cross(const AcReal3& a, const AcReal3& b)
{ {
@@ -63,9 +60,9 @@ vec_norm(const AcReal3& a)
AcReal3 c; AcReal3 c;
AcReal norm = dot(a, a); AcReal norm = dot(a, a);
c.x = a.x/sqrt(norm); c.x = a.x / sqrt(norm);
c.y = a.y/sqrt(norm); c.y = a.y / sqrt(norm);
c.z = a.z/sqrt(norm); c.z = a.z / sqrt(norm);
return c; return c;
} }
@@ -75,9 +72,9 @@ vec_multi_scal(const AcReal scal, const AcReal3& a)
{ {
AcReal3 c; AcReal3 c;
c.x = a.x*scal; c.x = a.x * scal;
c.y = a.y*scal; c.y = a.y * scal;
c.z = a.z*scal; c.z = a.z * scal;
return c; return c;
} }
@@ -86,95 +83,95 @@ vec_multi_scal(const AcReal scal, const AcReal3& a)
AcReal3 AcReal3
helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) helical_forcing_k_generator(const AcReal kmax, const AcReal kmin)
{ {
AcReal phi, theta, kk; //Spherical direction coordinates AcReal phi, theta, kk; // Spherical direction coordinates
AcReal3 k_force; //forcing wave vector AcReal3 k_force; // forcing wave vector
AcReal delta_k = kmax - kmin; AcReal delta_k = kmax - kmin;
// Generate vector in spherical coordinates // Generate vector in spherical coordinates
phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
theta = AcReal(M_PI)*get_random_number_01(); theta = AcReal(M_PI) * get_random_number_01();
kk = delta_k*get_random_number_01() + kmin; kk = delta_k * get_random_number_01() + kmin;
// Cast into Cartesian form // Cast into Cartesian form
k_force = (AcReal3){kk*sin(theta)*cos(phi), k_force = (AcReal3){kk * sin(theta) * cos(phi), //
kk*sin(theta)*sin(phi), kk * sin(theta) * sin(phi), //
kk*cos(theta) }; kk * cos(theta)};
//printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z); // printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
//Round the numbers. In that way k(x/y/z) will get complete waves. // Round the numbers. In that way k(x/y/z) will get complete waves.
k_force.x = round(k_force.x); k_force.y = round(k_force.y); k_force.z = round(k_force.z); k_force.x = round(k_force.x);
k_force.y = round(k_force.y);
k_force.z = round(k_force.z);
//printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z); // printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x,
// k_force.y, k_force.z);
return k_force; return k_force;
} }
//Generate the unit perpendicular unit vector e required for helical forcing // Generate the unit perpendicular unit vector e required for helical forcing
//Addapted from Pencil code forcing.f90 hel_vec() subroutine. // Addapted from Pencil code forcing.f90 hel_vec() subroutine.
void void
helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force)
{ {
AcReal3 k_cross_e = cross(k_force, *e_force); AcReal3 k_cross_e = cross(k_force, *e_force);
k_cross_e = vec_norm(k_cross_e); k_cross_e = vec_norm(k_cross_e);
AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e);
k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); k_cross_k_cross_e = vec_norm(k_cross_k_cross_e);
AcReal phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); AcReal phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); AcReal3 ee_tmp1 = vec_multi_scal(cos(phi), k_cross_e);
AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e);
*e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z};
ee_tmp1.y + ee_tmp2.y,
ee_tmp1.z + ee_tmp2.z};
} }
//PC Manual Eq. 223 // PC Manual Eq. 223
void void
helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force,
const AcReal3 e_force, const AcReal relhel) const AcReal3 e_force, const AcReal relhel)
{ {
// k dot e // k dot e
AcReal3 kdote; AcReal3 kdote;
kdote.x = k_force.x * e_force.x; kdote.x = k_force.x * e_force.x;
kdote.y = k_force.y * e_force.y; kdote.y = k_force.y * e_force.y;
kdote.z = k_force.z * e_force.z; kdote.z = k_force.z * e_force.z;
// k cross e // k cross e
AcReal3 k_cross_e; AcReal3 k_cross_e;
k_cross_e.x=k_force.y*e_force.z-k_force.z*e_force.y; k_cross_e.x = k_force.y * e_force.z - k_force.z * e_force.y;
k_cross_e.y=k_force.z*e_force.x-k_force.x*e_force.z; k_cross_e.y = k_force.z * e_force.x - k_force.x * e_force.z;
k_cross_e.z=k_force.x*e_force.y-k_force.y*e_force.x; k_cross_e.z = k_force.x * e_force.y - k_force.y * e_force.x;
// k cross k cross e // k cross k cross e
AcReal3 k_cross_k_cross_e; AcReal3 k_cross_k_cross_e;
k_cross_k_cross_e.x=k_force.y*k_cross_e.z-k_force.z*k_cross_e.y; k_cross_k_cross_e.x = k_force.y * k_cross_e.z - k_force.z * k_cross_e.y;
k_cross_k_cross_e.y=k_force.z*k_cross_e.x-k_force.x*k_cross_e.z; k_cross_k_cross_e.y = k_force.z * k_cross_e.x - k_force.x * k_cross_e.z;
k_cross_k_cross_e.z=k_force.x*k_cross_e.y-k_force.y*k_cross_e.x; k_cross_k_cross_e.z = k_force.x * k_cross_e.y - k_force.y * k_cross_e.x;
// abs(k) // abs(k)
AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); AcReal kabs = sqrt(k_force.x * k_force.x + k_force.y * k_force.y + k_force.z * k_force.z);
AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*kabs AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs *
*sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); sqrt(kabs * kabs -
(kdote.x * kdote.x + kdote.y * kdote.y + kdote.z * kdote.z));
//MV: I suspect there is a typo in the Pencil Code manual! // MV: I suspect there is a typo in the Pencil Code manual!
//*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator,
// -relhel*kabs*k_cross_e.y/denominator, // -relhel*kabs*k_cross_e.y/denominator,
// -relhel*kabs*k_cross_e.z/denominator}; // -relhel*kabs*k_cross_e.z/denominator};
//*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator,
// k_cross_k_cross_e.y/denominator, // k_cross_k_cross_e.y/denominator,
// k_cross_k_cross_e.z/denominator}; // k_cross_k_cross_e.z/denominator};
// See PC forcing.f90 forcing_hel_both() // See PC forcing.f90 forcing_hel_both()
*ff_hel_re = (AcReal3){kabs*k_cross_e.x/denominator, *ff_hel_re = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y,
kabs*k_cross_e.y, kabs * k_cross_e.z};
kabs*k_cross_e.z};
*ff_hel_im = (AcReal3){relhel*k_cross_k_cross_e.x/denominator, *ff_hel_im = (AcReal3){relhel * k_cross_k_cross_e.x / denominator, relhel * k_cross_k_cross_e.y,
relhel*k_cross_k_cross_e.y, relhel * k_cross_k_cross_e.z};
relhel*k_cross_k_cross_e.z};
} }