89 lines
2.1 KiB
Common Lisp
Raw Normal View History

2016-07-04 14:28:02 -07:00
/*
* H update equations
*
2016-08-03 14:51:25 -07:00
* Template parameters:
* mu False if (mu == 1) everywhere
* pmc False if no PMC anywhere
* common_cl Rendered code from common.cl
*
* Arguments:
* ctype *E E-field
* ctype *H H-field
* ctype *inv_mu 1/mu (at H-field locations)
* char *pmc Boolean mask denoting presence of PMC (at H-field locations)
* ctype *inv_dex 1/dx_e (complex cell widths for x direction at E locations)
* ctype *inv_dey 1/dy_e (complex cell widths for y direction at E locations)
* ctype *inv_dez 1/dz_e (complex cell widths for z direction at E locations)
*
2016-07-04 14:28:02 -07:00
*/
2016-08-03 14:51:25 -07:00
{{common_cl}}
2016-07-04 14:28:02 -07:00
////////////////////////////////////////////////////////////////////////////
2016-08-03 14:51:25 -07:00
__global ctype *inv_mu_x = inv_mu + XX;
__global ctype *inv_mu_y = inv_mu + YY;
__global ctype *inv_mu_z = inv_mu + ZZ;
2016-07-04 14:28:02 -07:00
2016-08-03 15:03:47 -07:00
__global char *pmc_x = pmc + XX;
__global char *pmc_y = pmc + YY;
__global char *pmc_z = pmc + ZZ;
2016-07-04 14:28:02 -07:00
2016-08-03 14:51:25 -07:00
//Update H components; set them to 0 if PMC is enabled at that location.
2017-05-20 21:30:32 -07:00
//Mu division and PMC conditional are only included if {mu} and {pmc} are true
2016-07-04 14:28:02 -07:00
{% if pmc -%}
2016-08-03 14:51:25 -07:00
if (pmc_x[i] != 0) {
Hx[i] = zero;
2016-07-04 14:28:02 -07:00
} else
{%- endif -%}
{
ctype Dzy = mul(sub(Ez[i + py], Ez[i]), inv_dey[y]);
ctype Dyz = mul(sub(Ey[i + pz], Ey[i]), inv_dez[z]);
2016-08-03 14:51:25 -07:00
ctype x_curl = sub(Dzy, Dyz);
2016-07-04 14:28:02 -07:00
2017-05-20 21:30:32 -07:00
{%- if mu %}
2016-08-03 14:51:25 -07:00
Hx[i] = mul(inv_mu_x[i], x_curl);
2017-05-20 21:30:32 -07:00
{%- else %}
2016-08-03 14:51:25 -07:00
Hx[i] = x_curl;
2016-07-04 14:28:02 -07:00
{%- endif %}
}
{% if pmc -%}
2016-08-03 14:51:25 -07:00
if (pmc_y[i] != 0) {
Hy[i] = zero;
2016-07-04 14:28:02 -07:00
} else
{%- endif -%}
{
ctype Dxz = mul(sub(Ex[i + pz], Ex[i]), inv_dez[z]);
ctype Dzx = mul(sub(Ez[i + px], Ez[i]), inv_dex[x]);
2016-08-03 14:51:25 -07:00
ctype y_curl = sub(Dxz, Dzx);
2016-07-04 14:28:02 -07:00
2017-05-20 21:30:32 -07:00
{%- if mu %}
2016-08-03 14:51:25 -07:00
Hy[i] = mul(inv_mu_y[i], y_curl);
2017-05-20 21:30:32 -07:00
{%- else %}
2016-08-03 14:51:25 -07:00
Hy[i] = y_curl;
2016-07-04 14:28:02 -07:00
{%- endif %}
}
{% if pmc -%}
2016-08-03 14:51:25 -07:00
if (pmc_z[i] != 0) {
Hz[i] = zero;
2016-07-04 14:28:02 -07:00
} else
{%- endif -%}
{
ctype Dyx = mul(sub(Ey[i + px], Ey[i]), inv_dex[x]);
ctype Dxy = mul(sub(Ex[i + py], Ex[i]), inv_dey[y]);
2016-08-03 14:51:25 -07:00
ctype z_curl = sub(Dyx, Dxy);
2016-07-04 14:28:02 -07:00
2017-05-20 21:30:32 -07:00
{%- if mu %}
2016-08-03 14:51:25 -07:00
Hz[i] = mul(inv_mu_z[i], z_curl);
2017-05-20 21:30:32 -07:00
{%- else %}
2016-08-03 14:51:25 -07:00
Hz[i] = z_curl;
2016-07-04 14:28:02 -07:00
{%- endif %}
}
/*
* End H update equations
2016-07-04 19:30:27 -07:00
*/