Special field updates for r = 0 and m = 0, ±1 in cylindrical coordinates #2538
Conversation
oskooi
changed the title
|
Added a new function As shown below, this fix produces the correct results whereby the DFT fields for
The caveat is there is still one case which does show a discrepancy:
|
I suspect that this is an unrelated problem with interpolation. |
src/meep_internals.hpp
Outdated
| #define STEP_UPDATE_EDHB0(f, fc, gv, is, ie, g, u, s, chi2, chi3, fw, dsigw, sigw, kapw) \ | ||
| do { \ | ||
| step_update_EDHB0(f, fc, gv, is, ie, g, u, s, chi2, chi3, fw, dsigw, sigw, kapw); \ | ||
| } while (0) |
There was a problem hiding this comment.
Why do we need a macro here rather than just calling step_update_EDHB0 directly?
There was a problem hiding this comment.
Actually, this line
Line 37 in e4504ed
is generating a
stride1 version called (confusingly) STEP_UPDATE_EDHB_stride10, which we should call from the macro similar to STEP_UPDATE_EDHB
|
I assume it still fixes #2182? |
The |
|
Perhaps there is separate bug in the step-curl update for Lines 313 to 347 in e4504ed |
|
Would be good to see Ep (before and after this PR) for m=1. Without this PR, it seems like Ep at r=0 would never be updated, which seems like it would screw up the Hr update at r=0, so I'm confused about why removing the Ep update would help. (Unless Ep is getting updated somewhere else?) |
|
with 2538
before 2538
|
|
This is confusing to me because it doesn't seem like Ep should be updated at all at r=0 without this patch. (The previous update_EDHB line skips r=0 for Ep because of the little_owned_corner0.) Maybe put in some printf calls to see where the Ep is being changed for r=0? |
|
Another quick check is to try setting μ = 1 + 1e-14 (i.e. force it to store and update H separately), in case the missing H update at r=0 (which affects the r=0 PML) is actually the real source of #2182. |
Looks like the reason the DFT and N2F fields for I inspected the time-domain (and not the DFT) fields at Summary of Results
Notes
diff --git a/src/update_eh.cpp b/src/update_eh.cpp
index 31e1c9e6..9b6d0c88 100644
--- a/src/update_eh.cpp
+++ b/src/update_eh.cpp
@@ -207,6 +207,20 @@ bool fields_chunk::update_eh(field_type ft, bool skip_w_components) {
}
}
}
+
+ if (ec == Ep || ec == Hr) {
+ realnum *E = f[ec][cmp];
+ const int yee_idx = gv.yee_index(ec);
+ const int sR = gv.stride(R), nZ = gv.num_direction(Z);
+ bool zero_everywhere = true;
+ for (int iZ = 0; iZ < nZ; iZ++) {
+ const int i = yee_idx + iZ - sR;
+ if (E[i] != 0)
+ zero_everywhere = false;
+ }
+ master_printf("zero-at-r0?, %s, %d, %d\n", component_name(ec), cmp, zero_everywhere);
+ }
+
}
}
}It seems therefore that this PR is producing the expected results (in the non-PML region anyway). However, we still need to explain why the |
I checked and found the ownership isn't changed by #1895, and #2182 was present before #1895. |
|
A couple of changes:
This PR does fix a bug for the There are two options for how to proceed: (1) we merge this PR as-is and then work on a bug fix for |
|
As discussed today, we went over the r=0 updates at Lines 313 to 346 in 955dcaa and so far they look correct, but we didn't check the PML terms yet. Would be good to compare those to the ones (for r≠0) in step_generic.cpp, since the PML equations should be the same (they don't involve spatial derivatives). |
|
I also verified that we are not adding the Lines 156 to 157 in eb8ff3a Lines 281 to 282 in eb8ff3a In the first instance, the loop over the Line 174 in eb8ff3a Line 182 in eb8ff3a |
The special step-curl update equation with Specifically, in Lines 334 to 341 in eb8ff3a The key line number is 339 involving updating For comparison, in Lines 197 to 204 in eb8ff3a Only the second term matches (The step-curl equation for |
|
I have tried to add the correct step-db update equations for However, there seems to be a bug somewhere in this implementation since the fields blow up for the unit test for |
|
|
There is a bug in the special step-db update equations at Clearly, the results for the cylindrical coordinates is wrong (the extraction efficiency cannot be larger than one). I went over the update equations again and compared them to the "most general case" of Lines 218 to 230 in 95ddacd |
|
There was a bug in the sign of the curl terms and strides which involved a simple fix: |
|
I think the old If a test is failing now and was passing with the old code in |
|
It seems there was a bug in the |
oskooi
changed the title
… Poynting flux of nearfields to avoid DFT interpolation to voxel center
|
In the failing unit test for the extraction efficiency which involves computing the radiated flux using I think the change from We can probably go ahead and merge this PR now and then file a bug report that we can address with a separate PR. |
|
It looks like that the raw Dp fields at r=0 actually do increase with resolution. At t=5, the fields are In addition, the fields oscillates twice as fast. For the same period of time=5, I have 500 snapshots of fields at resolution=25 and 1000 snapshots of fields at resolution=50. Therefore, frames 498-500 at resolution 25 should correspond to frames 996-1000 at resolution 50. However, plotting those frames, it looks like that the fields are oscillating at the rate of time stepping: 
|
|
@oskooi pointed out that, when resolution is doubled, pixel size is halved, so time-domain fields getting doubled would make sense. The evolution of fields in time is reasonable and consistent with what we expect from a correct time stepping. The next step is to check that there is no jump in the raw field as r approaches 0. |
For reference, due to the finite discretization a Lines 480 to 484 in dad3e52 Doubling the grid resolution should therefore result in a doubling of the field amplitudes in the time domain. This is what we are observing here. |
I realized that I was saving both real and imaginary parts of the fields together. If I focus on the real part (i.e. frames with odd number index), the fields evolution would make sense. Here is an animation for the real part of Dp at r=0 for t=0 to 5: r0_anim.mp4 |
|
On the other hand, using |
…r source at r > 0
|
Instead of a bug in the In the failing unit test of the extraction efficiency, if the position of the source is shifted by a voxel away from # Because (1) Er is not defined at r=0 on the Yee grid, and (2) there
# seems to be a bug in the restriction of an Er point source at r = 0,
# the source is placed at r = ~Δr (just outside the first voxel).
# This incurs a small error which decreases linearly with resolution.
src_pt = mp.Vector3(1.5 / self.resolution, 0, -0.5 * sz + h * dmat)Note that this hack only seems to work if the source position in |
|
Might be worth trying an integrated source, since that adds the current in a completely different place. Might help to determine whether this is problem with the Er field update. |
|
One thing to make sure of is that we are doing the transformation across r=0 correctly when interpolating the source in src_vol_chunkloop. Two things to try:
|
|
If the problem is with sources that touch the r=0 pixel, would be good to check whether the problem is simply an overall normalization or something deeper. One way to check this is to look at the spectrum, e.g. radiated power, for some nontrivial structure (e.g. a waveguide or sphere), and compare it for a source at r=0 and a source at r=1.5Δx. They should be the same, up to a first-order O(Δx) error, but if there is an overall normalization problem at r=0 they should still be nearly proportional to one another. |
|
Plot of the ratio of the radiated flux in air for the LED structure for two configurations of
|
|
This looks just like a scale factor (i.e. the ratio is nearly constant over most of the bandwidth). Presumably as you increase the resolution it should get more and more flat. That means that it's not a problem with the timestepping or anything like that — it's just an issue with the normalization when the source is placed. |
|
Here is a plot of Er field after 100 wavelength with source at r=0:
And here is a plot with source at r=1.5/Δr:
|
|
Some related problem might occur for the So, for example, if you did the extraction-efficiency tutorial but normalized by a |
I tried placing the This finding is unexpected because based on our discussions last week we had thought the bug in the @mochen, would you update your field plots above for an |
Here are the plots of fields with colorbars. I fixed their ranges to be the same: Er source at r=0: Er source at r=1.5Δr: |
I added the The output for the For both cases, there are only 2 grid points at ± 1 pixel in Δz. There is no restriction of the
The bug seems to affect only
|
|
1.5Δr will be exactly on the Yee grid point for Er, so you expect only 2 current points. It would be good to try 2.0Δr, in which case we expect 4 points. It would also be useful to add some print statements to
(A quick hack to avoid losing the r<0 contributions for an Er source at r=0 would be to double the source amplitude. But we should really fix the |
Findings:
In these four cases (output shown below), the sum of the amplitudes of the grid points is always 625.00. This indicates "charge conservation." These results suggest that the restriction for However, the sum of the amplitudes for an
|
|
Additional output from items 1-3 listed above. The main finding is that
For all five test cases, the output is always: |
|
Results from a separate experiment involving an The output is organized into four columns: Resolutions tested were 25, 50, 100, and 200 pixels/μm.
Cases 6-8 produce the correct result for the extraction efficiency because they do not involve restricting to |
|
Results for the experiment above but computing the total flux using The ratio of the radiated to the total flux is converging with resolution to the same constant value for all dipole positions.
|
|
Since the remaining problems are simply an overall normalization factor for Er delta function sources at (or near) r=0, I think we can go ahead an merge. |
Fixes #2489.
While #2383 fixed a bug for the$z$ -PML at $r = 0$ for $m = \pm 1$ , it seems to have introduced a separate bug for the $m = 0$ case. This is because #2383 actually never tested the $m = 0$ case which is why it went undetected until #2489.
This PR reverts the change to$m = 0$ . It is not yet clear to me why this fixes the issue reported in #2489 but the test results suggests it does.
update_eh.cppin #2383 but only forA new unit test for$m = 0$ (currently missing) based on the results in #2489 will also be added to this PR in a subsequent commit.
Note: two test cases in$m = 0$ required adjusting the tolerance of their error check.
test_adjoint_cyl.pyinvolving