From fad0af143a8640ef534f4c43872842ffb9634289 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Mon, 27 Nov 2023 21:29:55 -0800 Subject: [PATCH 1/2] two improvements to BFAST API --- python/examples/refl_angular_bfast.ipynb | 6 ++---- python/simulation.py | 9 +++------ python/tests/test_refl_angular.py | 12 +++++------- src/fields.cpp | 17 +++++++---------- src/meep.hpp | 12 +++++------- src/step_db.cpp | 7 +++++-- 6 files changed, 27 insertions(+), 36 deletions(-) diff --git a/python/examples/refl_angular_bfast.ipynb b/python/examples/refl_angular_bfast.ipynb index 483f6df17..b41b2e7a7 100644 --- a/python/examples/refl_angular_bfast.ipynb +++ b/python/examples/refl_angular_bfast.ipynb @@ -80,8 +80,7 @@ " k_point=mp.Vector3(),\n", " dimensions=dimensions,\n", " resolution=resolution,\n", - " need_bfast = True,\n", - " bfast_k_bar = (np.sin(theta_r),0,0), #sets the angle of the incident wave\n", + " bfast_scaled_k = (np.sin(theta_r),0,0), #sets the angle of the incident wave\n", " Courant = Courant \n", " )\n", "\n", @@ -116,8 +115,7 @@ " k_point=mp.Vector3(),\n", " dimensions=dimensions,\n", " resolution=resolution,\n", - " need_bfast = True,\n", - " bfast_k_bar = (np.sin(theta_r),0,0),\n", + " bfast_scaled_k = (np.sin(theta_r),0,0),\n", " Courant = Courant \n", " )\n", "\n", diff --git a/python/simulation.py b/python/simulation.py index 1abd3c995..91ce950de 100644 --- a/python/simulation.py +++ b/python/simulation.py @@ -1230,8 +1230,7 @@ def __init__( force_complex_fields: bool = False, default_material: Medium = mp.Medium(), m: float = 0, - need_bfast: bool = False, - bfast_k_bar: Vector3Type = (0, 0, 0), + bfast_scaled_k: Optional[Vector3Type] = None, k_point: Union[Vector3Type, bool] = False, kz_2d: str = "complex", extra_materials: Optional[List[Medium]] = None, @@ -1529,8 +1528,7 @@ def __init__( self.last_eps_filename = "" self.output_h5_hook = lambda fname: False self.interactive = False - self.need_bfast = need_bfast - self.bfast_k_bar = bfast_k_bar + self.bfast_scaled_k = (0, 0, 0) if bfast_scaled_k is None else bfast_scaled_k self.is_cylindrical = False self.material_function = material_function self.epsilon_func = epsilon_func @@ -2479,8 +2477,7 @@ def init_sim(self): not self.accurate_fields_near_cylorigin, self.loop_tile_base_db, self.loop_tile_base_eh, - self.need_bfast, - self.bfast_k_bar, + self.bfast_scaled_k, ) if self.force_all_components and self.dimensions != 1: diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index eb86f953e..d796cc57b 100644 --- a/python/tests/test_refl_angular.py +++ b/python/tests/test_refl_angular.py @@ -47,13 +47,13 @@ def reflectance_angular( theta_rad = math.radians(theta_deg) if need_bfast: - bfast_k_bar = (self.n1 * np.sin(theta_rad), 0, 0) + bfast_scaled_k = (self.n1 * np.sin(theta_rad), 0, 0) - Courant = (1 - bfast_k_bar[0]) / 3**0.5 + Courant = (1 - bfast_scaled_k[0]) / 3**0.5 k = mp.Vector3() else: - bfast_k_bar = (0, 0, 0) + bfast_scaled_k = (0, 0, 0) Courant = 0.5 @@ -89,8 +89,7 @@ def reflectance_angular( sources=sources, boundary_layers=pml_layers, k_point=k, - need_bfast=need_bfast, - bfast_k_bar=bfast_k_bar, + bfast_scaled_k=bfast_scaled_k, Courant=Courant, ) @@ -126,8 +125,7 @@ def reflectance_angular( sources=sources, boundary_layers=pml_layers, k_point=k, - need_bfast=need_bfast, - bfast_k_bar=bfast_k_bar, + bfast_scaled_k=bfast_scaled_k, Courant=Courant, geometry=geometry, ) diff --git a/src/fields.cpp b/src/fields.cpp index 93d79fa17..4d76c5026 100644 --- a/src/fields.cpp +++ b/src/fields.cpp @@ -30,11 +30,10 @@ using namespace std; namespace meep { fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylorigin, - int loop_tile_base_db, int loop_tile_base_eh, bool need_bfast, - std::vector bfast_k_bar) + int loop_tile_base_db, int loop_tile_base_eh, std::vector bfast_scaled_k) : S(s->S), gv(s->gv), user_volume(s->user_volume), v(s->v), m(m), beta(beta), loop_tile_base_db(loop_tile_base_db), loop_tile_base_eh(loop_tile_base_eh), - working_on(×_spent), need_bfast(need_bfast), bfast_k_bar(bfast_k_bar) { + working_on(×_spent), bfast_scaled_k(bfast_scaled_k) { shared_chunks = s->shared_chunks; components_allocated = false; synchronized_magnetic_fields = 0; @@ -61,7 +60,7 @@ fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylori chunks = new fields_chunk_ptr[num_chunks]; for (int i = 0; i < num_chunks; i++) chunks[i] = new fields_chunk(s->chunks[i], outdir, m, beta, zero_fields_near_cylorigin, i, - loop_tile_base_db, need_bfast, bfast_k_bar); + loop_tile_base_db, bfast_scaled_k); FOR_FIELD_TYPES(ft) { typedef realnum *realnum_ptr; comm_blocks[ft] = new realnum_ptr[num_chunks * num_chunks]; @@ -94,8 +93,7 @@ fields::fields(const fields &thef) outdir = new char[strlen(thef.outdir) + 1]; strcpy(outdir, thef.outdir); m = thef.m; - need_bfast = thef.need_bfast; - bfast_k_bar = thef.bfast_k_bar; + bfast_scaled_k = thef.bfast_scaled_k; beta = thef.beta; phasein_time = thef.phasein_time; for (int d = 0; d < 5; d++) { @@ -243,9 +241,9 @@ void check_tiles(grid_volume gv, const std::vector &gvs) { fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, double beta, bool zero_fields_near_cylorigin, int chunkidx, int loop_tile_base_db, - bool need_bfast, std::vector bfast_k_bar) + std::vector bfast_scaled_k) : gv(the_s->gv), v(the_s->v), m(m), zero_fields_near_cylorigin(zero_fields_near_cylorigin), - beta(beta), need_bfast(need_bfast), bfast_k_bar(bfast_k_bar) { + beta(beta), bfast_scaled_k(bfast_scaled_k) { s = the_s; chunk_idx = chunkidx; s->refcount++; @@ -307,8 +305,7 @@ fields_chunk::fields_chunk(const fields_chunk &thef, int chunkidx) : gv(thef.gv) s->refcount++; outdir = thef.outdir; m = thef.m; - need_bfast = thef.need_bfast; - bfast_k_bar = thef.bfast_k_bar; + bfast_scaled_k = thef.bfast_scaled_k; zero_fields_near_cylorigin = thef.zero_fields_near_cylorigin; beta = thef.beta; new_s = thef.new_s; diff --git a/src/meep.hpp b/src/meep.hpp index 8fc10fc5b..371de5a13 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1502,8 +1502,7 @@ class fields_chunk { volume v; double m; // angular dependence in cyl. coords bool zero_fields_near_cylorigin; // fields=0 m pixels near r=0 for stability - bool need_bfast; - std::vector bfast_k_bar; + std::vector bfast_scaled_k; double beta; int is_real; std::vector sources[NUM_FIELD_TYPES]; @@ -1514,7 +1513,7 @@ class fields_chunk { fields_chunk(structure_chunk *, const char *outdir, double m, double beta, bool zero_fields_near_cylorigin, int chunkidx, int loop_tile_base_db, - bool need_bfast, std::vector bfast_k_bar); + std::vector bfast_scaled_k); fields_chunk(const fields_chunk &, int chunkidx); ~fields_chunk(); @@ -1742,8 +1741,7 @@ class fields { grid_volume gv, user_volume; volume v; double m; - bool need_bfast; - std::vector bfast_k_bar; + std::vector bfast_scaled_k; double beta; int t, phasein_time, is_real; std::complex k[5], eikna[5]; @@ -1755,8 +1753,8 @@ class fields { // fields.cpp methods: fields(structure *, double m = 0, double beta = 0, bool zero_fields_near_cylorigin = true, - int loop_tile_base_db = 0, int loop_tile_base_eh = 0, bool need_bfast = false, - std::vector bfast_k_bar = {0, 0, 0}); + int loop_tile_base_db = 0, int loop_tile_base_eh = 0, + std::vector bfast_scaled_k = {0, 0, 0}); fields(const fields &); ~fields(); bool equal_layout(const fields &f) const; diff --git a/src/step_db.cpp b/src/step_db.cpp index 21b227ffb..24c27b62c 100644 --- a/src/step_db.cpp +++ b/src/step_db.cpp @@ -62,6 +62,7 @@ bool fields_chunk::step_db(field_type ft) { realnum *f_p = have_p ? f[c_p][cmp] : NULL; realnum *f_m = have_m ? f[c_m][cmp] : NULL; realnum *the_f = f[cc][cmp]; + bool need_bfast = (bfast_scaled_k[0] || bfast_scaled_k[1] || bfast_scaled_k[2]) ? 1 : 0; if (dsig != NO_DIRECTION && s->conductivity[cc][d_c] && !f_cond[cc][cmp]) { f_cond[cc][cmp] = new realnum[gv.ntot()]; @@ -126,8 +127,10 @@ bool fields_chunk::step_db(field_type ft) { s->conductivity[cc][d_c], s->condinv[cc][d_c], f_cond[cc][cmp]); if (need_bfast) { - realnum k1 = have_m ? bfast_k_bar[component_index(c_m)] : 0; // puts k1 in direction of g2 - realnum k2 = have_p ? bfast_k_bar[component_index(c_p)] : 0; // puts k2 in direction of g1 + realnum k1 = + have_m ? bfast_scaled_k[component_index(c_m)] : 0; // puts k1 in direction of g2 + realnum k2 = + have_p ? bfast_scaled_k[component_index(c_p)] : 0; // puts k2 in direction of g1 if (ft == D_stuff) { k1 = -k1; k2 = -k2; From ee9e68d0cffe1f97609162aefa18ddfe09bdd122 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Tue, 28 Nov 2023 23:43:45 +0000 Subject: [PATCH 2/2] rename need_bfast to use_bfast and remove ternary operator --- python/tests/test_refl_angular.py | 12 ++++++------ src/step_db.cpp | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py index d796cc57b..dd3221b44 100644 --- a/python/tests/test_refl_angular.py +++ b/python/tests/test_refl_angular.py @@ -30,13 +30,13 @@ def setUpClass(cls): cls.num_freq = 11 def reflectance_angular( - self, theta_deg: float, need_bfast: bool + self, theta_deg: float, use_bfast: bool ) -> Tuple[List, List, np.ndarray]: """Computes properties of the incident and reflected planewave. Args: theta_deg: angle of incident planewave. - need_bfast: whether to use the same angle for the incident planewave + use_bfast: whether to use the same angle for the incident planewave for all frequencies. If False, the incident angle is frequency dependent. @@ -46,7 +46,7 @@ def reflectance_angular( """ theta_rad = math.radians(theta_deg) - if need_bfast: + if use_bfast: bfast_scaled_k = (self.n1 * np.sin(theta_rad), 0, 0) Courant = (1 - bfast_scaled_k[0]) / 3**0.5 @@ -142,7 +142,7 @@ def reflectance_angular( reflectance = -np.array(flux_monitor_flux) / np.array(empty_flux) - if need_bfast: + if use_bfast: theta_in_rad = [theta_rad] * self.num_freq else: # Returns the angle of the incident planewave in medium n1 based @@ -154,12 +154,12 @@ def reflectance_angular( return freqs, theta_in_rad, reflectance @parameterized.parameterized.expand([(0, False), (20.6, False), (35.7, True)]) - def test_reflectance_angular(self, theta_deg: float, need_bfast: bool): + def test_reflectance_angular(self, theta_deg: float, use_bfast: bool): ( frequency_meep, theta_in_rad_meep, reflectance_meep, - ) = self.reflectance_angular(theta_deg, need_bfast) + ) = self.reflectance_angular(theta_deg, use_bfast) # Returns angle of refracted planewave in medium n2 given # an incident planewave in medium n1 at angle theta_in_rad. diff --git a/src/step_db.cpp b/src/step_db.cpp index 24c27b62c..67e4dbebb 100644 --- a/src/step_db.cpp +++ b/src/step_db.cpp @@ -62,7 +62,7 @@ bool fields_chunk::step_db(field_type ft) { realnum *f_p = have_p ? f[c_p][cmp] : NULL; realnum *f_m = have_m ? f[c_m][cmp] : NULL; realnum *the_f = f[cc][cmp]; - bool need_bfast = (bfast_scaled_k[0] || bfast_scaled_k[1] || bfast_scaled_k[2]) ? 1 : 0; + bool use_bfast = bfast_scaled_k[0] || bfast_scaled_k[1] || bfast_scaled_k[2]; if (dsig != NO_DIRECTION && s->conductivity[cc][d_c] && !f_cond[cc][cmp]) { f_cond[cc][cmp] = new realnum[gv.ntot()]; @@ -73,7 +73,7 @@ bool fields_chunk::step_db(field_type ft) { memcpy(f_u[cc][cmp], the_f, gv.ntot() * sizeof(realnum)); allocated_u = true; } - if (need_bfast && !f_bfast[cc][cmp]) { + if (use_bfast && !f_bfast[cc][cmp]) { f_bfast[cc][cmp] = new realnum[gv.ntot()]; memset(f_bfast[cc][cmp], 0, sizeof(realnum) * gv.ntot()); } @@ -126,7 +126,7 @@ bool fields_chunk::step_db(field_type ft) { 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]); - if (need_bfast) { + if (use_bfast) { realnum k1 = have_m ? bfast_scaled_k[component_index(c_m)] : 0; // puts k1 in direction of g2 realnum k2 =