From 4e48e241549af56d1a107ea735f10edd207a8035 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 20 Apr 2023 08:09:45 -0700 Subject: [PATCH 1/2] place cylindrical dipole sources with r weight --- src/meep.hpp | 7 +++++++ src/sources.cpp | 29 +++++++++++++++++++++++++++-- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index a1f728d40..0c5e96dde 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1524,6 +1524,10 @@ class fields_chunk { } double last_source_time(); + + // sources.cpp + std::complex sum_sources(field_type ft); + // monitor.cpp std::complex get_field(component, const ivec &) const; @@ -1761,6 +1765,9 @@ class fields { void change_m(double new_m); bool has_nonlinearities(bool parallel = true) const; + // sources.coo + std::complex sum_sources(field_type ft); + // time.cpp std::vector time_spent_on(time_sink sink); double mean_time_spent_on(time_sink); diff --git a/src/sources.cpp b/src/sources.cpp index e4ef6d1bd..93a8b80d0 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -270,7 +270,12 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is loc += shift * (0.5 * inva); vec rel_loc = loc - data->center; - amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc); + + if (fc->gv.dim == Dcyl) + amps_array[idx_vol] = + IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2) * amp * data->A(rel_loc); + else + amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc); // check for invalid sources at r=0 in cylindrical coordinates if (fc->gv.dim == Dcyl && loc.r() == 0 && amps_array[idx_vol] != 0.0) { @@ -479,7 +484,7 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w data.A = A ? A : one; data.amp = amp; LOOP_OVER_DIRECTIONS(gv.dim, d) { - if (where.in_direction(d) == 0.0 && !nosize_direction(d)) // delta-fun + if (gv.dim != Dcyl && where.in_direction(d) == 0.0 && !nosize_direction(d)) // delta-fun data.amp *= gv.a; // correct units for J delta-function amplitude } sources = src.add_to(sources, &data.src); @@ -755,4 +760,24 @@ diffractedplanewave::diffractedplanewave(int g_[3], double axis_[3], std::comple p = p_; }; +std::complex fields::sum_sources(field_type ft) { + std::complex total_sources = 0.0; + + for (int i = 0; i < num_chunks; i++) + if (chunks[i]->is_mine()) total_sources += chunks[i]->sum_sources(ft); + + return sum_to_all(total_sources); +} + +std::complex fields_chunk::sum_sources(field_type ft) { + std::complex total_sources = 0.0; + for (const src_vol &sv : sources[ft]) { + for (size_t j = 0; j < sv.num_points(); j++) { + meep::master_printf("src %f\n", sv.amplitude_at(j)); + total_sources += sv.amplitude_at(j); + } + } + return total_sources; +} + } // namespace meep From 825ee616e6738aa9333c1fd4231d3bb627c12dc1 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 20 Apr 2023 11:08:20 -0700 Subject: [PATCH 2/2] compensate for r --- src/sources.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/sources.cpp b/src/sources.cpp index 93a8b80d0..f0c257d3d 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -271,6 +271,7 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is vec rel_loc = loc - data->center; + // in cylindrical coordinates, we need to account for the r term in the integral. if (fc->gv.dim == Dcyl) amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2) * amp * data->A(rel_loc); @@ -483,9 +484,16 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w src_vol_chunkloop_data data; data.A = A ? A : one; data.amp = amp; + // correct units for J delta-function amplitude LOOP_OVER_DIRECTIONS(gv.dim, d) { - if (gv.dim != Dcyl && where.in_direction(d) == 0.0 && !nosize_direction(d)) // delta-fun - data.amp *= gv.a; // correct units for J delta-function amplitude + if (where.in_direction(d) == 0.0 && !nosize_direction(d)) { // delta-fun + if (d == R && where.center().in_direction(d) > 0) + // in cylindrical coordinates, we need to scale by $r$ for a dipole. + data.amp *= 1.0 / where.center().in_direction(d); + else if (gv.dim != Dcyl) + // we handle the resolution scaling at the chunk level for cylindrical coordinates. + data.amp *= gv.a; + } } sources = src.add_to(sources, &data.src); data.center = (where.get_min_corner() + where.get_max_corner()) * 0.5; @@ -760,6 +768,11 @@ diffractedplanewave::diffractedplanewave(int g_[3], double axis_[3], std::comple p = p_; }; +/* +Occasionally, it's useful to check the sum of all the actual source amplitudes +(e.g. to ensure that the boundary weights are being computed properly). +`sum_sources()` provides a convenient interface to do just that. +*/ std::complex fields::sum_sources(field_type ft) { std::complex total_sources = 0.0;