Adapted the old example of helical forcing with profiles to conform with the revised syntax

This commit is contained in:
jpekkila
2019-10-07 19:43:25 +03:00
parent 0cc5bdaa08
commit ea4438f331

View File

@@ -0,0 +1,55 @@
// PC-style helical forcing with profiles
uniform Vector AC_kk;
uniform Vector AC_coef1;
uniform Vector AC_coef2;
uniform Vector AC_coef3;
uniform Vector AC_fda;
uniform Scalar AC_phase;
uniform Scalar AC_fact;
uniform Scalar AC_k1_ff;
uniform ScalarArray AC_profx;
uniform ScalarArray AC_profy;
uniform ScalarArray AC_profz;
uniform ScalarArray AC_profx_hel;
uniform ScalarArray AC_profy_hel;
uniform ScalarArray AC_profz_hel;
uniform int AC_iforcing_zsym;
// PC-style helical forcing with support for profiles
Vector
pcforcing(int3 vertexIdx, int3 globalVertexIdx, Scalar dt, ScalarArray profx, ScalarArray profy,
ScalarArray profz, ScalarArray profx_hel, ScalarArray profy_hel, ScalarArray profz_hel)
{
Vector pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz};
Complex fx = AC_fact * exp(Complex(0, AC_kk.x * AC_k1_ff * pos.z + AC_phase));
Complex fy = exp(Complex(0, AC_kk.y * AC_k1_ff * pos.y));
Complex fz;
if (AC_iforcing_zsym == 0) {
fz = exp(Complex(0.0, AC_kk.z * AC_k1_ff * pos.z));
}
else if (AC_iforcing_zsym == 1) {
fz = Complex(cos(AC_kk.z * AC_k1_ff * pos.z), 0);
}
else if (AC_iforcing_zsym == -1) {
fz = Complex(sin(AC_kk.z * AC_k1_ff * pos.z), 0);
}
else {
// Failure
}
Complex fxyz = fx * fy * fz;
// TODO recheck indices
Scalar force_ampl = profx[vertexIdx.x - NGHOST] * profy[vertexIdx.y] * profz[vertexIdx.z];
Scalar prof_hel_ampl = profx_hel[vertexIdx.x - NGHOST] * profy_hel[vertexIdx.y] *
profz_hel[vertexIdx.z];
return force_ampl * AC_fda * (Complex(AC_coef1.x, prof_hel_ampl * AC_coef2.x) * fxyz).x;
}