diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas index 53b5a45..5c4f14a 100644 --- a/acc/mhd_solver/stencil_assembly.sas +++ b/acc/mhd_solver/stencil_assembly.sas @@ -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; diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 813f51a..ee8656e 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -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 diff --git a/acc/src/acc.l b/acc/src/acc.l index e68fe8b..b180ae8 100644 --- a/acc/src/acc.l +++ b/acc/src/acc.l @@ -15,6 +15,8 @@ L [a-zA-Z_] "void" { return VOID; } /* Rest of the types inherited from C */ "int" { return INT; } "int3" { return INT3; } +"ScalarField" { return SCALAR; } +"VectorField" { return VECTOR; } "Kernel" { return KERNEL; } /* Function specifiers */ "Preprocessed" { return PREPROCESSED; } diff --git a/acc/src/acc.y b/acc/src/acc.y index db49225..7138b1b 100644 --- a/acc/src/acc.y +++ b/acc/src/acc.y @@ -20,7 +20,7 @@ int yyget_lineno(); %token VOID INT INT3 %token IF ELSE FOR WHILE ELIF %token LEQU LAND LOR LLEQU -%token KERNEL PREPROCESSED +%token KERNEL PREPROCESSED %token INPLACE_INC INPLACE_DEC %% @@ -72,11 +72,11 @@ selection_statement: IF expression else_selection_statement ; else_selection_statement: compound_statement { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); } - | compound_statement elif_selection_statement { $$ = astnode_create(NODE_UNKNOWN, $1, $2); } + | compound_statement elif_selection_statement { $$ = astnode_create(NODE_UNKNOWN, $1, $2); } | compound_statement ELSE compound_statement { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = ELSE; } ; -elif_selection_statement: ELIF expression else_selection_statement { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = ELIF; } +elif_selection_statement: ELIF expression else_selection_statement { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = ELIF; } ; iteration_statement: WHILE expression compound_statement { $$ = astnode_create(NODE_UNKNOWN, $2, $3); $$->prefix = WHILE; } @@ -101,8 +101,9 @@ exec_statement: declaration ; assignment: declaration '=' expression { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '='; } + | declaration '(' expression_list ')' { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '('; $$->postfix = ')'; } // C++ style initializer | expression '=' expression { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '='; } - ; + ; return_statement: /* Empty */ { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); } | expression { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); } @@ -126,7 +127,7 @@ array_declaration: identifier '[' ']' | identifier '[' expression ']' { $$ = astnode_create(NODE_UNKNOWN, $1, $3); $$->infix = '['; $$->postfix = ']'; } ; -type_declaration: type_specifier { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); } +type_declaration: type_specifier { $$ = astnode_create(NODE_UNKNOWN, $1, NULL); } | type_qualifier type_specifier { $$ = astnode_create(NODE_UNKNOWN, $1, $2); } ; @@ -181,7 +182,7 @@ binary_operator: '+' | '/' { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; } | '*' { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; } | '<' { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; } - | '>' { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; } + | '>' { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); $$->infix = yytext[0]; } | LEQU { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); } | LAND { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); } | LOR { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); } diff --git a/acc/src/code_generator.c b/acc/src/code_generator.c index 93ca42f..648fc9b 100644 --- a/acc/src/code_generator.c +++ b/acc/src/code_generator.c @@ -229,8 +229,11 @@ translate_latest_symbol(void) // IN / OUT else if (symbol->type != SYMBOLTYPE_FUNCTION_PARAMETER && (symbol->type_qualifier == IN || symbol->type_qualifier == OUT)) { - const char* inout_type_qualifier = "static __device__ const auto"; - printf("%s %s%s", inout_type_qualifier, inout_name_prefix, symbol_table[handle].identifier); + + printf("static __device__ const %s %s%s", symbol->type_specifier == SCALAR ? "int" : "int3", + inout_name_prefix, symbol_table[handle].identifier); + if (symbol->type_specifier == VECTOR) + printf(" = make_int3"); } // OTHER else { diff --git a/src/core/kernels/integration.cuh b/src/core/kernels/integration.cuh index fec3b95..50db9ae 100644 --- a/src/core/kernels/integration.cuh +++ b/src/core/kernels/integration.cuh @@ -46,6 +46,9 @@ IDX(const int3 idx) return DEVICE_VTXBUF_IDX(idx.x, idx.y, idx.z); } +#define make_int3(a, b, c) \ + (int3) { (int)a, (int)b, (int)c } + static __forceinline__ AcMatrix create_rotz(const AcReal radians) {