diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc index ecac783..3c4c714 100644 --- a/src/standalone/model/host_forcing.cc +++ b/src/standalone/model/host_forcing.cc @@ -28,17 +28,14 @@ #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 get_random_number_01() { - //TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ - return AcReal(rand())/AcReal(RAND_MAX); + // TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ + return AcReal(rand()) / AcReal(RAND_MAX); } - - AcReal3 cross(const AcReal3& a, const AcReal3& b) { @@ -63,9 +60,9 @@ vec_norm(const AcReal3& a) AcReal3 c; AcReal norm = dot(a, a); - c.x = a.x/sqrt(norm); - c.y = a.y/sqrt(norm); - c.z = a.z/sqrt(norm); + c.x = a.x / sqrt(norm); + c.y = a.y / sqrt(norm); + c.z = a.z / sqrt(norm); return c; } @@ -75,9 +72,9 @@ vec_multi_scal(const AcReal scal, const AcReal3& a) { AcReal3 c; - c.x = a.x*scal; - c.y = a.y*scal; - c.z = a.z*scal; + c.x = a.x * scal; + c.y = a.y * scal; + c.z = a.z * scal; return c; } @@ -86,95 +83,95 @@ vec_multi_scal(const AcReal scal, const AcReal3& a) AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) { - AcReal phi, theta, kk; //Spherical direction coordinates - AcReal3 k_force; //forcing wave vector + AcReal phi, theta, kk; // Spherical direction coordinates + AcReal3 k_force; // forcing wave vector AcReal delta_k = kmax - kmin; // Generate vector in spherical coordinates - phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); - theta = AcReal(M_PI)*get_random_number_01(); - kk = delta_k*get_random_number_01() + kmin; + phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01(); + theta = AcReal(M_PI) * get_random_number_01(); + kk = delta_k * get_random_number_01() + kmin; // Cast into Cartesian form - k_force = (AcReal3){kk*sin(theta)*cos(phi), - kk*sin(theta)*sin(phi), - kk*cos(theta) }; + k_force = (AcReal3){kk * sin(theta) * cos(phi), // + kk * sin(theta) * sin(phi), // + 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. - k_force.x = round(k_force.x); k_force.y = round(k_force.y); k_force.z = round(k_force.z); + // 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); - //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; } -//Generate the unit perpendicular unit vector e required for helical forcing -//Addapted from Pencil code forcing.f90 hel_vec() subroutine. +// Generate the unit perpendicular unit vector e required for helical forcing +// Addapted from Pencil code forcing.f90 hel_vec() subroutine. void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) { - AcReal3 k_cross_e = cross(k_force, *e_force); - k_cross_e = vec_norm(k_cross_e); + AcReal3 k_cross_e = cross(k_force, *e_force); + k_cross_e = vec_norm(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); - 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_tmp2 = vec_multi_scal(sin(phi), 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(); + AcReal3 ee_tmp1 = vec_multi_scal(cos(phi), 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, - ee_tmp1.y + ee_tmp2.y, - ee_tmp1.z + ee_tmp2.z}; + *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z}; } -//PC Manual Eq. 223 +// PC Manual Eq. 223 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) { - // k dot e + // k dot e AcReal3 kdote; kdote.x = k_force.x * e_force.x; kdote.y = k_force.y * e_force.y; kdote.z = k_force.z * e_force.z; // 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.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; + AcReal3 k_cross_e; + 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.z = k_force.x * e_force.y - k_force.y * e_force.x; // 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.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; + 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.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; // 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 - *sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); + AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs * + 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! - //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, + // MV: I suspect there is a typo in the Pencil Code manual! + //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, // -relhel*kabs*k_cross_e.y/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.z/denominator}; // See PC forcing.f90 forcing_hel_both() - *ff_hel_re = (AcReal3){kabs*k_cross_e.x/denominator, - kabs*k_cross_e.y, - kabs*k_cross_e.z}; + *ff_hel_re = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y, + kabs * k_cross_e.z}; - *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.z}; + *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.z}; }