diff --git a/src/meep_internals.hpp b/src/meep_internals.hpp index 051b5e0d3..0bb56b5a5 100644 --- a/src/meep_internals.hpp +++ b/src/meep_internals.hpp @@ -95,18 +95,19 @@ symmetry r_to_minus_r_symmetry(int m); void step_curl(realnum *f, component c, const realnum *g1, const realnum *g2, ptrdiff_t s1, ptrdiff_t s2, // strides for g1/g2 shift - const grid_volume &gv, realnum dtdx, direction dsig, const realnum *sig, + const grid_volume &gv, const ivec is, const ivec ie, + realnum dtdx, direction dsig, const realnum *sig, const realnum *kap, const realnum *siginv, realnum *fu, direction dsigu, const realnum *sigu, const realnum *kapu, const realnum *siginvu, realnum dt, const realnum *cnd, const realnum *cndinv, realnum *fcnd); -void step_update_EDHB(realnum *f, component fc, const grid_volume &gv, const realnum *g, +void step_update_EDHB(realnum *f, component fc, const grid_volume &gv, const ivec is, const ivec ie, const realnum *g, const realnum *g1, const realnum *g2, const realnum *u, const realnum *u1, const realnum *u2, ptrdiff_t s, ptrdiff_t s1, ptrdiff_t s2, const realnum *chi2, const realnum *chi3, realnum *fw, direction dsigw, const realnum *sigw, const realnum *kapw); -void step_beta(realnum *f, component c, const realnum *g, const grid_volume &gv, realnum betadt, +void step_beta(realnum *f, component c, const realnum *g, const grid_volume &gv, const ivec is, const ivec ie, realnum betadt, direction dsig, const realnum *siginv, realnum *fu, direction dsigu, const realnum *siginvu, const realnum *cndinv, realnum *fcnd); @@ -114,18 +115,19 @@ void step_beta(realnum *f, component c, const realnum *g, const grid_volume &gv, void step_curl_stride1(realnum *f, component c, const realnum *g1, const realnum *g2, ptrdiff_t s1, ptrdiff_t s2, // strides for g1/g2 shift - const grid_volume &gv, realnum dtdx, direction dsig, const realnum *sig, + const grid_volume &gv, const ivec is, const ivec ie, + realnum dtdx, direction dsig, const realnum *sig, const realnum *kap, const realnum *siginv, realnum *fu, direction dsigu, const realnum *sigu, const realnum *kapu, const realnum *siginvu, realnum dt, const realnum *cnd, const realnum *cndinv, realnum *fcnd); -void step_update_EDHB_stride1(realnum *f, component fc, const grid_volume &gv, const realnum *g, +void step_update_EDHB_stride1(realnum *f, component fc, const grid_volume &gv, const ivec is, const ivec ie, const realnum *g, const realnum *g1, const realnum *g2, const realnum *u, const realnum *u1, const realnum *u2, ptrdiff_t s, ptrdiff_t s1, ptrdiff_t s2, const realnum *chi2, const realnum *chi3, realnum *fw, direction dsigw, const realnum *sigw, const realnum *kapw); -void step_beta_stride1(realnum *f, component c, const realnum *g, const grid_volume &gv, +void step_beta_stride1(realnum *f, component c, const realnum *g, const grid_volume &gv, const ivec is, const ivec ie, realnum betadt, direction dsig, const realnum *siginv, realnum *fu, direction dsigu, const realnum *siginvu, const realnum *cndinv, realnum *fcnd); @@ -135,34 +137,34 @@ void step_beta_stride1(realnum *f, component c, const realnum *g, const grid_vol which allow gcc (and possibly other compilers) to do additional optimizations, especially loop vectorization */ -#define STEP_CURL(f, c, g1, g2, s1, s2, gv, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, kapu, \ +#define STEP_CURL(f, c, g1, g2, s1, s2, gv, is, ie, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, kapu, \ siginvu, dt, cnd, cndinv, fcnd) \ do { \ if (LOOPS_ARE_STRIDE1(gv)) \ - step_curl_stride1(f, c, g1, g2, s1, s2, gv, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, \ + step_curl_stride1(f, c, g1, g2, s1, s2, gv, is, ie, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, \ kapu, siginvu, dt, cnd, cndinv, fcnd); \ else \ - step_curl(f, c, g1, g2, s1, s2, gv, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, kapu, \ + step_curl(f, c, g1, g2, s1, s2, gv, is, ie, dtdx, dsig, sig, kap, siginv, fu, dsigu, sigu, kapu, \ siginvu, dt, cnd, cndinv, fcnd); \ } while (0) -#define STEP_UPDATE_EDHB(f, fc, gv, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, sigw, \ +#define STEP_UPDATE_EDHB(f, fc, gv, is, ie, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, sigw, \ kapw) \ do { \ if (LOOPS_ARE_STRIDE1(gv)) \ - step_update_EDHB_stride1(f, fc, gv, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, \ + step_update_EDHB_stride1(f, fc, gv, is, ie, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, \ sigw, kapw); \ else \ - step_update_EDHB(f, fc, gv, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, sigw, \ + step_update_EDHB(f, fc, gv, is, ie, g, g1, g2, u, u1, u2, s, s1, s2, chi2, chi3, fw, dsigw, sigw, \ kapw); \ } while (0) -#define STEP_BETA(f, c, g, gv, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd) \ +#define STEP_BETA(f, c, g, gv, is, ie, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd) \ do { \ if (LOOPS_ARE_STRIDE1(gv)) \ - step_beta_stride1(f, c, g, gv, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd); \ + step_beta_stride1(f, c, g, gv, is, ie, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd); \ else \ - step_beta(f, c, g, gv, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd); \ + step_beta(f, c, g, gv, is, ie, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd); \ } while (0) // analytical Green's functions from near2far.cpp, which we might want to expose someday diff --git a/src/step_db.cpp b/src/step_db.cpp index 378cc5d2e..5a31f93bc 100644 --- a/src/step_db.cpp +++ b/src/step_db.cpp @@ -115,7 +115,9 @@ bool fields_chunk::step_db(field_type ft) { default: meep::abort("bug - non-cylindrical field component in Dcyl"); } - STEP_CURL(the_f, cc, f_p, f_m, stride_p, stride_m, gv, Courant, dsig, s->sig[dsig], + STEP_CURL(the_f, cc, f_p, f_m, stride_p, stride_m, + gv, gv.little_owned_corner0(cc), gv.big_corner(), + Courant, dsig, s->sig[dsig], s->kap[dsig], s->siginv[dsig], f_u[cc][cmp], dsigu, s->sig[dsigu], s->kap[dsigu], s->siginv[dsigu], dt, s->conductivity[cc][d_c], s->condinv[cc][d_c], f_cond[cc][cmp]); @@ -146,7 +148,8 @@ bool fields_chunk::step_db(field_type ft) { const direction dsigu = s->sigsize[dsigu0] > 1 ? dsigu0 : NO_DIRECTION; const realnum betadt = 2 * pi * beta * dt * (d_c == X ? +1 : -1) * (f[c_g][1 - cmp] ? (ft == D_stuff ? -1 : +1) * (2 * cmp - 1) : 1); - STEP_BETA(the_f, cc, g, gv, betadt, dsig, s->siginv[dsig], f_u[cc][cmp], dsigu, + STEP_BETA(the_f, cc, g, gv, gv.little_owned_corner0(cc), gv.big_corner(), + betadt, dsig, s->siginv[dsig], f_u[cc][cmp], dsigu, s->siginv[dsigu], s->condinv[cc][d_c], f_cond[cc][cmp]); } diff --git a/src/step_generic.cpp b/src/step_generic.cpp index 62e06be57..66e79a745 100644 --- a/src/step_generic.cpp +++ b/src/step_generic.cpp @@ -65,7 +65,8 @@ namespace meep { */ void step_curl(RPR f, component c, const RPR g1, const RPR g2, ptrdiff_t s1, ptrdiff_t s2, // strides for g1/g2 shift - const grid_volume &gv, realnum dtdx, direction dsig, const RPR sig, const RPR kap, + const grid_volume &gv, const ivec is, const ivec ie, + realnum dtdx, direction dsig, const RPR sig, const RPR kap, const RPR siginv, RPR fu, direction dsigu, const RPR sigu, const RPR kapu, const RPR siginvu, realnum dt, const RPR cnd, const RPR cndinv, RPR fcnd) { if (!g1) { // swap g1 and g2 @@ -86,25 +87,25 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, if (cnd) { realnum dt2 = dt * 0.5; if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] = ((1 - dt2 * cnd[i]) * f[i] - dtdx * (g1[i + s1] - g1[i] + g2[i] - g2[i + s2])) * cndinv[i]; } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] = ((1 - dt2 * cnd[i]) * f[i] - dtdx * (g1[i + s1] - g1[i])) * cndinv[i]; } } } else { // no conductivity if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] -= dtdx * (g1[i + s1] - g1[i] + g2[i] - g2[i + s2]); } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { f[i] -= dtdx * (g1[i + s1] - g1[i]); } + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] -= dtdx * (g1[i + s1] - g1[i]); } } } } @@ -113,7 +114,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, if (cnd) { realnum dt2 = dt * 0.5; if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum fprev = fu[i]; fu[i] = @@ -123,7 +124,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum fprev = fu[i]; fu[i] = ((1 - dt2 * cnd[i]) * fprev - dtdx * (g1[i + s1] - g1[i])) * cndinv[i]; @@ -133,7 +134,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } else { // no conductivity if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum fprev = fu[i]; fu[i] -= dtdx * (g1[i + s1] - g1[i] + g2[i] - g2[i + s2]); @@ -141,7 +142,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum fprev = fu[i]; fu[i] -= dtdx * (g1[i + s1] - g1[i]); @@ -157,7 +158,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, if (cnd) { realnum dt2 = dt * 0.5; if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; realnum fcnd_prev = fcnd[i]; fcnd[i] = @@ -167,7 +168,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; realnum fcnd_prev = fcnd[i]; fcnd[i] = ((1 - dt2 * cnd[i]) * fcnd[i] - dtdx * (g1[i + s1] - g1[i])) * cndinv[i]; @@ -177,14 +178,14 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } else { // no conductivity (other than PML conductivity) if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; f[i] = ((kap[k] - sig[k]) * f[i] - dtdx * (g1[i + s1] - g1[i] + g2[i] - g2[i + s2])) * siginv[k]; } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; f[i] = ((kap[k] - sig[k]) * f[i] - dtdx * (g1[i + s1] - g1[i])) * siginv[k]; } @@ -197,7 +198,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, realnum dt2 = dt * 0.5; if (g2) { //////////////////// MOST GENERAL CASE ////////////////////// - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum fprev = fu[i]; @@ -211,7 +212,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, ///////////////////////////////////////////////////////////// } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum fprev = fu[i]; @@ -224,7 +225,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } else { // no conductivity (other than PML conductivity) if (g2) { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum fprev = fu[i]; @@ -234,7 +235,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, } } else { - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum fprev = fu[i]; @@ -251,7 +252,7 @@ void step_curl(RPR f, component c, const RPR g1, const RPR g2, and/or PML). This is used in 2d calculations to add an exp(i beta z) time dependence, which gives an additional i \beta \hat{z} \times cross-product in the curl equations. */ -void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum betadt, +void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, const ivec is, const ivec ie, realnum betadt, direction dsig, const RPR siginv, RPR fu, direction dsigu, const RPR siginvu, const RPR cndinv, RPR fcnd) { if (!g) return; @@ -261,7 +262,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b KSTRIDE_DEF(dsigu, ku, gv.little_owned_corner0(c)); if (cndinv) { // conductivity + PML //////////////////// MOST GENERAL CASE ////////////////////// - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum df; @@ -273,7 +274,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b ///////////////////////////////////////////////////////////// } else { // PML only - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; DEF_ku; realnum df; @@ -284,7 +285,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b } else { // PML in f, no fu if (cndinv) { // conductivity + PML - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; realnum dfcnd = betadt * g[i] * cndinv[i]; fcnd[i] += dfcnd; @@ -292,7 +293,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b } } else { // PML only - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_k; f[i] += betadt * g[i] * siginv[k]; } @@ -303,7 +304,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b if (dsigu != NO_DIRECTION) { // fu, no PML in f KSTRIDE_DEF(dsigu, ku, gv.little_owned_corner0(c)); if (cndinv) { // conductivity, no PML - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum df; fu[i] += (df = betadt * g[i] * cndinv[i]); @@ -311,7 +312,7 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b } } else { // no conductivity or PML - LOOP_OVER_VOL_OWNED0(gv, c, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_ku; realnum df; fu[i] += (df = betadt * g[i]); @@ -321,10 +322,10 @@ void step_beta(RPR f, component c, const RPR g, const grid_volume &gv, realnum b } else { // no PML, no fu if (cndinv) { // conductivity, no PML - LOOP_OVER_VOL_OWNED0(gv, c, i) { f[i] += betadt * g[i] * cndinv[i]; } + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] += betadt * g[i] * cndinv[i]; } } else { // no conductivity or PML - LOOP_OVER_VOL_OWNED0(gv, c, i) { f[i] += betadt * g[i]; } + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] += betadt * g[i]; } } } } @@ -361,7 +362,7 @@ inline realnum calc_nonlinear_u(const realnum Dsqr, const realnum Di, const real */ -void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, const RPR g1, +void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const ivec is, const ivec ie, const RPR g, const RPR g1, const RPR g2, const RPR u, const RPR u1, const RPR u2, ptrdiff_t s, ptrdiff_t s1, ptrdiff_t s2, const RPR chi2, const RPR chi3, RPR fw, direction dsigw, const RPR sigw, const RPR kapw) { @@ -385,7 +386,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c if (u1 && u2) { // 3x3 off-diagonal u if (chi3) { //////////////////// MOST GENERAL CASE ////////////////////// - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum g2s = g2[i] + g2[i + s] + g2[i - s2] + g2[i + (s - s2)]; realnum gs = g[i]; @@ -400,7 +401,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c ///////////////////////////////////////////////////////////// } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; DEF_kw; @@ -412,7 +413,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } else if (u1) { // 2x2 off-diagonal u if (chi3) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum gs = g[i]; realnum us = u[i]; @@ -424,7 +425,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; DEF_kw; @@ -440,7 +441,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c else { // diagonal u if (chi3) { if (g1 && g2) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum g2s = g2[i] + g2[i + s] + g2[i - s2] + g2[i + (s - s2)]; realnum gs = g[i]; @@ -453,7 +454,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else if (g1) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum gs = g[i]; realnum us = u[i]; @@ -468,7 +469,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c meep::abort("bug - didn't swap off-diagonal terms!?"); } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; DEF_kw; @@ -479,7 +480,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else if (u) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; DEF_kw; @@ -489,7 +490,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { DEF_kw; realnum fwprev = fw[i], kapwkw = kapw[kw], sigwkw = sigw[kw]; fw[i] = g[i]; @@ -501,7 +502,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c else { /////////////// no PML (no fw) /////////////////// if (u1 && u2) { // 3x3 off-diagonal u if (chi3) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum g2s = g2[i] + g2[i + s] + g2[i - s2] + g2[i + (s - s2)]; realnum gs = g[i]; @@ -512,7 +513,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; f[i] = (gs * us + OFFDIAG(u1, g1, s1) + OFFDIAG(u2, g2, s2)); @@ -521,7 +522,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } else if (u1) { // 2x2 off-diagonal u if (chi3) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum gs = g[i]; realnum us = u[i]; @@ -530,7 +531,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; f[i] = (gs * us + OFFDIAG(u1, g1, s1)); @@ -543,7 +544,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c else { // diagonal u if (chi3) { if (g1 && g2) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum g2s = g2[i] + g2[i + s] + g2[i - s2] + g2[i + (s - s2)]; realnum gs = g[i]; @@ -553,7 +554,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else if (g1) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum g1s = g1[i] + g1[i + s] + g1[i - s1] + g1[i + (s - s1)]; realnum gs = g[i]; realnum us = u[i]; @@ -565,7 +566,7 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c meep::abort("bug - didn't swap off-diagonal terms!?"); } else { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; f[i] = (gs * us) * calc_nonlinear_u(gs * gs, gs, us, chi2[i], chi3[i]); @@ -573,14 +574,14 @@ void step_update_EDHB(RPR f, component fc, const grid_volume &gv, const RPR g, c } } else if (u) { - LOOP_OVER_VOL_OWNED(gv, fc, i) { + LOOP_OVER_IVECS(gv, is, ie, i) { realnum gs = g[i]; realnum us = u[i]; f[i] = (gs * us); } } else - LOOP_OVER_VOL_OWNED(gv, fc, i) { f[i] = g[i]; } + LOOP_OVER_IVECS(gv, is, ie, i) { f[i] = g[i]; } } } } diff --git a/src/update_eh.cpp b/src/update_eh.cpp index 858e5f267..fd86ac993 100644 --- a/src/update_eh.cpp +++ b/src/update_eh.cpp @@ -165,7 +165,8 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) { } if (f[ec][cmp] != f[dc][cmp]) - STEP_UPDATE_EDHB(f[ec][cmp], ec, gv, dmp[dc][cmp], dmp[dc_1][cmp], dmp[dc_2][cmp], + STEP_UPDATE_EDHB(f[ec][cmp], ec, gv, gv.little_owned_corner(ec), gv.big_corner(), + dmp[dc][cmp], dmp[dc_1][cmp], dmp[dc_2][cmp], s->chi1inv[ec][d_ec], dmp[dc_1][cmp] ? s->chi1inv[ec][d_1] : NULL, dmp[dc_2][cmp] ? s->chi1inv[ec][d_2] : NULL, s_ec, s_1, s_2, s->chi2[ec], s->chi3[ec], f_w[ec][cmp], dsigw, s->sig[dsigw], s->kap[dsigw]);