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

Moved the old mhd solver to mhd_solver_DEPRECATED and replaced it with the new stencil_kernel.ac file

Browse Source
This commit is contained in:
jpekkila
2019-10-07 17:36:23 +03:00
parent 8c1e603a98
commit 6560be7056
4 changed files with 680 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

75
acc/mhd_solver/stencil_assembly.sas
View File

@@ -1,75 +0,0 @@
#include "stencil_definition.sdh"
Preprocessed Scalar
value(in ScalarField vertex)
{
return vertex[vertexIdx];
}
Preprocessed Vector
gradient(in ScalarField vertex)
{
return (Vector){derx(vertexIdx, vertex), dery(vertexIdx, vertex), derz(vertexIdx, vertex)};
}
#if LUPWD
Preprocessed Scalar
der6x_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsx;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 1, vertexIdx.y, vertexIdx.z]) -
Scalar(6.0) * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 2, vertexIdx.y, vertexIdx.z]) +
vertex[vertexIdx.x + 3, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 3, vertexIdx.y, vertexIdx.z])};
}
Preprocessed Scalar
der6y_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsy;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 1, vertexIdx.z]) -
Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 2, vertexIdx.z]) +
vertex[vertexIdx.x, vertexIdx.y + 3, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 3, vertexIdx.z])};
}
Preprocessed Scalar
der6z_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsz;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 1]) -
Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 2] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 2]) +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])};
}
#endif
Preprocessed Matrix
hessian(in ScalarField vertex)
{
Matrix hessian;
hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex),
derxz(vertexIdx, vertex)};
hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)};
hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(vertexIdx, vertex)};
return hessian;
}

137
acc/mhd_solver/stencil_definition.sdh
View File

@@ -1,137 +0,0 @@
#define LDENSITY (1)
#define LHYDRO (1)
#define LMAGNETIC (1)
#define LENTROPY (1)
#define LTEMPERATURE (0)
#define LFORCING (1)
#define LUPWD (1)
#define LSINK (0)
#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
// Int params
uniform int AC_max_steps;
uniform int AC_save_steps;
uniform int AC_bin_steps;
uniform int AC_bc_type;
uniform int AC_start_step;
// Real params
uniform Scalar AC_dt;
uniform Scalar AC_max_time;
// Spacing
uniform Scalar AC_dsx;
uniform Scalar AC_dsy;
uniform Scalar AC_dsz;
uniform Scalar AC_dsmin;
// physical grid
uniform Scalar AC_xlen;
uniform Scalar AC_ylen;
uniform Scalar AC_zlen;
uniform Scalar AC_xorig;
uniform Scalar AC_yorig;
uniform Scalar AC_zorig;
// Physical units
uniform Scalar AC_unit_density;
uniform Scalar AC_unit_velocity;
uniform Scalar AC_unit_length;
// properties of gravitating star
uniform Scalar AC_star_pos_x;
uniform Scalar AC_star_pos_y;
uniform Scalar AC_star_pos_z;
uniform Scalar AC_M_star;
// properties of sink particle
uniform Scalar AC_sink_pos_x;
uniform Scalar AC_sink_pos_y;
uniform Scalar AC_sink_pos_z;
uniform Scalar AC_M_sink;
uniform Scalar AC_M_sink_init;
uniform Scalar AC_M_sink_Msun;
uniform Scalar AC_soft;
uniform Scalar AC_accretion_range;
uniform Scalar AC_switch_accretion;
// Run params
uniform Scalar AC_cdt;
uniform Scalar AC_cdtv;
uniform Scalar AC_cdts;
uniform Scalar AC_nu_visc;
uniform Scalar AC_cs_sound;
uniform Scalar AC_eta;
uniform Scalar AC_mu0;
uniform Scalar AC_cp_sound;
uniform Scalar AC_gamma;
uniform Scalar AC_cv_sound;
uniform Scalar AC_lnT0;
uniform Scalar AC_lnrho0;
uniform Scalar AC_zeta;
uniform Scalar AC_trans;
// Other
uniform Scalar AC_bin_save_t;
// Initial condition params
uniform Scalar AC_ampl_lnrho;
uniform Scalar AC_ampl_uu;
uniform Scalar AC_angl_uu;
uniform Scalar AC_lnrho_edge;
uniform Scalar AC_lnrho_out;
// Forcing parameters. User configured.
uniform Scalar AC_forcing_magnitude;
uniform Scalar AC_relhel;
uniform Scalar AC_kmin;
uniform Scalar AC_kmax;
// Forcing parameters. Set by the generator.
uniform Scalar AC_forcing_phase;
uniform Scalar AC_k_forcex;
uniform Scalar AC_k_forcey;
uniform Scalar AC_k_forcez;
uniform Scalar AC_kaver;
uniform Scalar AC_ff_hel_rex;
uniform Scalar AC_ff_hel_rey;
uniform Scalar AC_ff_hel_rez;
uniform Scalar AC_ff_hel_imx;
uniform Scalar AC_ff_hel_imy;
uniform Scalar AC_ff_hel_imz;
// Additional helper params // (deduced from other params do not set these directly!)
uniform Scalar AC_G_const;
uniform Scalar AC_GM_star;
uniform Scalar AC_unit_mass;
uniform Scalar AC_sq2GM_star;
uniform Scalar AC_cs2_sound;
uniform Scalar AC_inv_dsx;
uniform Scalar AC_inv_dsy;
uniform Scalar AC_inv_dsz;
/*
* =============================================================================
* User-defined vertex buffers
* =============================================================================
*/
#if LENTROPY
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
uniform ScalarField VTXBUF_AX;
uniform ScalarField VTXBUF_AY;
uniform ScalarField VTXBUF_AZ;
uniform ScalarField VTXBUF_ENTROPY;
#elif LMAGNETIC
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
uniform ScalarField VTXBUF_AX;
uniform ScalarField VTXBUF_AY;
uniform ScalarField VTXBUF_AZ;
#elif LHYDRO
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
#else
uniform ScalarField VTXBUF_LNRHO;
#endif
#if LSINK
uniform ScalarField VTXBUF_ACCRETION;
#endif

256
acc/mhd_solver/stencil_process.sps → acc/mhd_solver/stencil_kernel.ac
View File

@@ -1,13 +1,190 @@
#include "stencil_definition.sdh"
#include <stdderiv.h>
Vector
value(in VectorField uu)
{
return (Vector){value(uu.x), value(uu.y), value(uu.z)};
}
#define LDENSITY (1)
#define LHYDRO (1)
#define LMAGNETIC (1)
#define LENTROPY (1)
#define LTEMPERATURE (0)
#define LFORCING (1)
#define LUPWD (1)
#define LSINK (0)
#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter
#define H_CONST (0) // TODO: make an actual config parameter
#define C_CONST (0) // TODO: make an actual config parameter
// Int params
uniform int AC_max_steps;
uniform int AC_save_steps;
uniform int AC_bin_steps;
uniform int AC_bc_type;
uniform int AC_start_step;
// Real params
uniform Scalar AC_dt;
uniform Scalar AC_max_time;
// Spacing
uniform Scalar AC_dsmin;
// physical grid
uniform Scalar AC_xlen;
uniform Scalar AC_ylen;
uniform Scalar AC_zlen;
uniform Scalar AC_xorig;
uniform Scalar AC_yorig;
uniform Scalar AC_zorig;
// Physical units
uniform Scalar AC_unit_density;
uniform Scalar AC_unit_velocity;
uniform Scalar AC_unit_length;
// properties of gravitating star
uniform Scalar AC_star_pos_x;
uniform Scalar AC_star_pos_y;
uniform Scalar AC_star_pos_z;
uniform Scalar AC_M_star;
// properties of sink particle
uniform Scalar AC_sink_pos_x;
uniform Scalar AC_sink_pos_y;
uniform Scalar AC_sink_pos_z;
uniform Scalar AC_M_sink;
uniform Scalar AC_M_sink_init;
uniform Scalar AC_M_sink_Msun;
uniform Scalar AC_soft;
uniform Scalar AC_accretion_range;
uniform Scalar AC_switch_accretion;
// Run params
uniform Scalar AC_cdt;
uniform Scalar AC_cdtv;
uniform Scalar AC_cdts;
uniform Scalar AC_nu_visc;
uniform Scalar AC_cs_sound;
uniform Scalar AC_eta;
uniform Scalar AC_mu0;
uniform Scalar AC_cp_sound;
uniform Scalar AC_gamma;
uniform Scalar AC_cv_sound;
uniform Scalar AC_lnT0;
uniform Scalar AC_lnrho0;
uniform Scalar AC_zeta;
uniform Scalar AC_trans;
// Other
uniform Scalar AC_bin_save_t;
// Initial condition params
uniform Scalar AC_ampl_lnrho;
uniform Scalar AC_ampl_uu;
uniform Scalar AC_angl_uu;
uniform Scalar AC_lnrho_edge;
uniform Scalar AC_lnrho_out;
// Forcing parameters. User configured.
uniform Scalar AC_forcing_magnitude;
uniform Scalar AC_relhel;
uniform Scalar AC_kmin;
uniform Scalar AC_kmax;
// Forcing parameters. Set by the generator.
uniform Scalar AC_forcing_phase;
uniform Scalar AC_k_forcex;
uniform Scalar AC_k_forcey;
uniform Scalar AC_k_forcez;
uniform Scalar AC_kaver;
uniform Scalar AC_ff_hel_rex;
uniform Scalar AC_ff_hel_rey;
uniform Scalar AC_ff_hel_rez;
uniform Scalar AC_ff_hel_imx;
uniform Scalar AC_ff_hel_imy;
uniform Scalar AC_ff_hel_imz;
// Additional helper params // (deduced from other params do not set these directly!)
uniform Scalar AC_G_const;
uniform Scalar AC_GM_star;
uniform Scalar AC_unit_mass;
uniform Scalar AC_sq2GM_star;
uniform Scalar AC_cs2_sound;
/*
* =============================================================================
* User-defined vertex buffers
* =============================================================================
*/
#if LENTROPY
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
uniform ScalarField VTXBUF_AX;
uniform ScalarField VTXBUF_AY;
uniform ScalarField VTXBUF_AZ;
uniform ScalarField VTXBUF_ENTROPY;
#elif LMAGNETIC
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
uniform ScalarField VTXBUF_AX;
uniform ScalarField VTXBUF_AY;
uniform ScalarField VTXBUF_AZ;
#elif LHYDRO
uniform ScalarField VTXBUF_LNRHO;
uniform ScalarField VTXBUF_UUX;
uniform ScalarField VTXBUF_UUY;
uniform ScalarField VTXBUF_UUZ;
#else
uniform ScalarField VTXBUF_LNRHO;
#endif
#if LSINK
uniform ScalarField VTXBUF_ACCRETION;
#endif
#if LUPWD
Scalar
Preprocessed Scalar
der6x_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsx;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 1, vertexIdx.y, vertexIdx.z]) -
Scalar(6.0) * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 2, vertexIdx.y, vertexIdx.z]) +
vertex[vertexIdx.x + 3, vertexIdx.y, vertexIdx.z] +
vertex[vertexIdx.x - 3, vertexIdx.y, vertexIdx.z])};
}
Preprocessed Scalar
der6y_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsy;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 1, vertexIdx.z]) -
Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 2, vertexIdx.z]) +
vertex[vertexIdx.x, vertexIdx.y + 3, vertexIdx.z] +
vertex[vertexIdx.x, vertexIdx.y - 3, vertexIdx.z])};
}
Preprocessed Scalar
der6z_upwd(in ScalarField vertex)
{
Scalar inv_ds = AC_inv_dsz;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 1]) -
Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 2] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 2]) +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])};
}
#endif
#if LUPWD
Device Scalar
upwd_der6(in VectorField uu, in ScalarField lnrho)
{
Scalar uux = fabs(value(uu).x);
@@ -17,21 +194,21 @@ upwd_der6(in VectorField uu, in ScalarField lnrho)
}
#endif
Matrix
Device Matrix
gradients(in VectorField uu)
{
return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
}
#if LSINK
Vector
Device Vector
sink_gravity(int3 globalVertexIdx){
int accretion_switch = int(AC_switch_accretion);
if (accretion_switch == 1){
Vector force_gravity;
const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Scalar sink_mass = AC_M_sink;
const Vector sink_pos = (Vector){AC_sink_pos_x,
AC_sink_pos_y,
@@ -55,7 +232,7 @@ sink_gravity(int3 globalVertexIdx){
#if LSINK
// Give Truelove density
Scalar
Device Scalar
truelove_density(in ScalarField lnrho){
const Scalar rho = exp(value(lnrho));
const Scalar Jeans_length_squared = (M_PI * AC_cs2_sound) / (AC_G_const * rho);
@@ -68,11 +245,11 @@ truelove_density(in ScalarField lnrho){
}
// This controls accretion of density/mass to the sink particle.
Scalar
Device Scalar
sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Vector sink_pos = (Vector){AC_sink_pos_x,
AC_sink_pos_y,
AC_sink_pos_z};
@@ -108,11 +285,11 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
}
// This controls accretion of velocity to the sink particle.
Vector
Device Vector
sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
const Vector grid_pos = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz};
const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Vector sink_pos = (Vector){AC_sink_pos_x,
AC_sink_pos_y,
AC_sink_pos_z};
@@ -153,7 +330,7 @@ sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
#endif
Scalar
Device Scalar
continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt)
{
return -dot(value(uu), gradient(lnrho))
@@ -170,7 +347,7 @@ continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar
#if LENTROPY
Vector
Device Vector
momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa, Scalar dt)
{
const Matrix S = stress_tensor(uu);
@@ -204,7 +381,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala
return mom;
}
#elif LTEMPERATURE
Vector
Device Vector
momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField tt)
{
Vector mom;
@@ -230,7 +407,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala
return mom;
}
#else
Vector
Device Vector
momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt)
{
Vector mom;
@@ -261,7 +438,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar d
}
#endif
Vector
Device Vector
induction(in VectorField uu, in VectorField aa)
{
// Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla
@@ -279,16 +456,15 @@ induction(in VectorField uu, in VectorField aa)
}
#if LENTROPY
Scalar
Device Scalar
lnT(in ScalarField ss, in ScalarField lnrho)
{
const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound +
return AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound +
(AC_gamma - Scalar(1.0)) * (value(lnrho) - AC_lnrho0);
return lnT;
}
// Nabla dot (K nabla T) / (rho T)
Scalar
Device Scalar
heat_conduction(in ScalarField ss, in ScalarField lnrho)
{
const Scalar inv_AC_cp_sound = AcReal(1.0) / AC_cp_sound;
@@ -306,13 +482,13 @@ heat_conduction(in ScalarField ss, in ScalarField lnrho)
return AC_cp_sound * chi * (first_term + dot(second_term, third_term));
}
Scalar
Device Scalar
heating(const int i, const int j, const int k)
{
return 1;
}
Scalar
Device Scalar
entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa)
{
const Matrix S = stress_tensor(uu);
@@ -328,7 +504,7 @@ entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorFie
#endif
#if LTEMPERATURE
Scalar
Device Scalar
heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
{
const Matrix S = stress_tensor(uu);
@@ -341,7 +517,7 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
#endif
#if LFORCING
Vector
Device Vector
simple_vortex_forcing(Vector a, Vector b, Scalar magnitude){
int accretion_switch = AC_switch_accretion;
@@ -351,7 +527,7 @@ Vector
return (Vector){0,0,0};
}
}
Vector
Device Vector
simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude){
int accretion_switch = AC_switch_accretion;
if (accretion_switch == 0){
@@ -363,7 +539,7 @@ Vector
// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable
// helicity
Vector
Device Vector
helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
{
// JP: This looks wrong:
@@ -397,7 +573,7 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
return force;
}
Vector
Device Vector
forcing(int3 globalVertexIdx, Scalar dt)
{
int accretion_switch = AC_switch_accretion;
@@ -406,9 +582,9 @@ forcing(int3 globalVertexIdx, Scalar dt)
Vector a = Scalar(0.5) * (Vector){globalGridN.x * AC_dsx,
globalGridN.y * AC_dsy,
globalGridN.z * AC_dsz}; // source (origin)
Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index)
Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz}; // sink (current index)
const Scalar cs2 = AC_cs2_sound;
const Scalar cs = sqrt(cs2);