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

Made the DSL syntax less confusing: Input and output arrays are now ScalarField and VectorFields instead of scalars and vectors. C++ initializers are now also possible, removing the need to declare Fields as int or int3 which was very confusing, like "what, you assing an int value to a real, what the &^%@?"

Browse Source
This commit is contained in:
jpekkila
2019-08-08 21:07:36 +03:00
parent 5397495496
commit b53cabbc44
6 changed files with 54 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
acc/mhd_solver/stencil_assembly.sas
View File

@@ -1,11 +1,11 @@
Preprocessed Scalar
value(in Scalar vertex)
value(in ScalarField vertex)
{
return vertex[vertexIdx];
}
Preprocessed Vector
gradient(in Scalar vertex)
gradient(in ScalarField vertex)
{
return (Vector){derx(vertexIdx, vertex),
dery(vertexIdx, vertex),
@@ -15,54 +15,54 @@ gradient(in Scalar vertex)
#if LUPWD
Preprocessed Scalar
der6x_upwd(in Scalar vertex)
der6x_upwd(in ScalarField vertex)
{
Scalar inv_ds = DCONST_REAL(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]
+ 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]
- 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]
+ vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])};
}
Preprocessed Scalar
der6y_upwd(in Scalar vertex)
der6y_upwd(in ScalarField vertex)
{
Scalar inv_ds = DCONST_REAL(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]
+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]
-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]
+ vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])};
}
Preprocessed Scalar
der6z_upwd(in Scalar vertex)
der6z_upwd(in ScalarField vertex)
{
Scalar inv_ds = DCONST_REAL(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]
+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]
-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]
+ vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])};
}
#endif
Preprocessed Matrix
hessian(in Scalar vertex)
hessian(in ScalarField vertex)
{
Matrix hessian;

44
acc/mhd_solver/stencil_process.sps
View File

@@ -25,14 +25,14 @@ uniform int nz;
Vector
value(in Vector uu)
value(in VectorField uu)
{
return (Vector){value(uu.x), value(uu.y), value(uu.z)};
}
#if LUPWD
Scalar
upwd_der6(in Vector uu, in Scalar lnrho)
upwd_der6(in VectorField uu, in ScalarField lnrho)
{
Scalar uux = fabs(value(uu).x);
Scalar uuy = fabs(value(uu).y);
@@ -42,13 +42,13 @@ upwd_der6(in Vector uu, in Scalar lnrho)
#endif
Matrix
gradients(in Vector uu)
gradients(in VectorField uu)
{
return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)};
}
Scalar
continuity(in Vector uu, in Scalar lnrho) {
continuity(in VectorField uu, in ScalarField lnrho) {
return -dot(value(uu), gradient(lnrho))
#if LUPWD
//This is a corrective hyperdiffusion term for upwinding.
@@ -59,7 +59,7 @@ continuity(in Vector uu, in Scalar lnrho) {
#if LENTROPY
Vector
momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) {
const Matrix S = stress_tensor(uu);
const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - lnrho0));
const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
@@ -82,7 +82,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) {
}
#elif LTEMPERATURE
Vector
momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) {
Vector mom;
const Matrix S = stress_tensor(uu);
@@ -103,7 +103,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar tt) {
}
#else
Vector
momentum(in Vector uu, in Scalar lnrho) {
momentum(in VectorField uu, in ScalarField lnrho) {
Vector mom;
const Matrix S = stress_tensor(uu);
@@ -126,7 +126,7 @@ momentum(in Vector uu, in Scalar lnrho) {
Vector
induction(in Vector uu, in Vector aa) {
induction(in VectorField uu, in VectorField aa) {
// Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla
// x A)) in order to avoid taking the first derivative twice (did the math,
// yes this actually works. See pg.28 in arXiv:astro-ph/0109497)
@@ -144,7 +144,7 @@ induction(in Vector uu, in Vector aa) {
#if LENTROPY
Scalar
lnT( in Scalar ss, in Scalar lnrho) {
lnT( in ScalarField ss, in ScalarField lnrho) {
const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound +
(gamma - Scalar(1.)) * (value(lnrho) - lnrho0);
return lnT;
@@ -152,7 +152,7 @@ lnT( in Scalar ss, in Scalar lnrho) {
// Nabla dot (K nabla T) / (rho T)
Scalar
heat_conduction( in Scalar ss, in Scalar lnrho) {
heat_conduction( in ScalarField ss, in ScalarField lnrho) {
const Scalar inv_cp_sound = AcReal(1.) / cp_sound;
const Vector grad_ln_chi = - gradient(lnrho);
@@ -174,7 +174,7 @@ heating(const int i, const int j, const int k) {
}
Scalar
entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) {
const Matrix S = stress_tensor(uu);
const Scalar inv_pT = Scalar(1.) / (exp(value(lnrho)) * exp(lnT(ss, lnrho)));
const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
@@ -191,7 +191,7 @@ entropy(in Scalar ss, in Vector uu, in Scalar lnrho, in Vector aa) {
#if LTEMPERATURE
Scalar
heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
{
const Matrix S = stress_tensor(uu);
const Scalar heat_diffusivity_k = 0.0008; //8e-4;
@@ -290,26 +290,26 @@ forcing(int3 globalVertexIdx, Scalar dt)
// Declare input and output arrays using locations specified in the
// array enum in astaroth.h
in Scalar lnrho = VTXBUF_LNRHO;
out Scalar out_lnrho = VTXBUF_LNRHO;
in ScalarField lnrho(VTXBUF_LNRHO);
out ScalarField out_lnrho(VTXBUF_LNRHO);
in Vector uu = (int3) {VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ};
out Vector out_uu = (int3) {VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ};
in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ);
#if LMAGNETIC
in Vector aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
out Vector out_aa = (int3) {VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ};
in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ);
#endif
#if LENTROPY
in Scalar ss = VTXBUF_ENTROPY;
out Scalar out_ss = VTXBUF_ENTROPY;
in ScalarField ss(VTXBUF_ENTROPY);
out ScalarField out_ss(VTXBUF_ENTROPY);
#endif
#if LTEMPERATURE
in Scalar tt = VTXBUF_TEMPERATURE;
out Scalar out_tt = VTXBUF_TEMPERATURE;
in ScalarField tt(VTXBUF_TEMPERATURE);
out ScalarField out_tt(VTXBUF_TEMPERATURE);
#endif
Kernel void