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

Scetching the helical forcing.

Browse Source 
Not the idead form. Not yet tested.
This commit is contained in:
Miikka Vaisala
2019-06-25 19:04:53 +08:00
parent a574d6e4c3
commit 8191c47fa0
2 changed files with 40 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
acc/mhd_solver/stencil_assembly.sas
View File

@@ -1,5 +1,5 @@
#define LUPWD (0) #define LUPWD (1)
Preprocessed Scalar Preprocessed Scalar

44
acc/mhd_solver/stencil_process.sps
View File

@@ -2,8 +2,8 @@
#define LENTROPY (1) #define LENTROPY (1)
#define LTEMPERATURE (0) #define LTEMPERATURE (0)
#define LGRAVITY (0) #define LGRAVITY (0)
#define LFORCING (0) #define LFORCING (1)
#define LUPWD (0) #define LUPWD (1)
// Declare uniforms (i.e. device constants) // Declare uniforms (i.e. device constants)
@@ -203,8 +203,12 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
} }
#endif #endif
#if LFORCING #if LFORCING
Vector Vector
simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
{ {
@@ -217,23 +221,51 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
} }
// The Pencil Code hel_vec(), manual Eq. 222, inspired forcing function with adjustable helicity
Vector
helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
{
Scalar cos_phi = cos(phi);
Scalar sin_phi = sin(phi);
Scalar cos_k_dox_x = cos(dot(k_force, xx));
Scalar sin_k_dox_x = sin(dot(k_force, xx));
Scalar real_comp = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi;
Scalar imag_comp = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi;
Vector force = (Vector){ ff_re.x*real_comp - ff_im.x*imag_comp,
ff_re.y*real_comp - ff_im.y*imag_comp,
ff_re.z*real_comp - ff_im.z*imag_comp};
return force;
}
Vector Vector
forcing(int3 globalVertexIdx) forcing(int3 globalVertexIdx)
{ {
Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx, Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx,
globalGrid.n.y * dsy, globalGrid.n.y * dsy,
globalGrid.n.z * dsz}; // source (origin) globalGrid.n.z * dsz}; // source (origin)
Vector b = (Vector){(globalVertexIdx.x - nx_min) * dsx, Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx,
(globalVertexIdx.y - ny_min) * dsy, (globalVertexIdx.y - ny_min) * dsy,
(globalVertexIdx.z - nz_min) * dsz}; // sink (current index) (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
//Placeholders until determined properly
Scalar magnitude = 0.05; Scalar magnitude = 0.05;
Vector k_force = (Vector){2.0, 0.0, 0.0};
Vector ff_re = (Vector){0.0, 0.5, 0.0};
Vector ff_im = (Vector){0.0, 0.8666, 0.0};
Scalar phase = Scalar(0.79);
//Determine that forcing funtion type at this point. //Determine that forcing funtion type at this point.
Vector c = simple_vortex_forcing(a, b, magnitude); //Vector force = simple_vortex_forcing(a, xx, magnitude);
//Vector c = simple_outward_flow_forcing(a, b, magnitude); //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phi);
if (is_valid(c)) { return c; } if (is_valid(force)) { return force; }
else { return (Vector){0, 0, 0}; } else { return (Vector){0, 0, 0}; }
} }
#endif // LFORCING #endif // LFORCING