bfast algorithm for fixed angle broadband

[Dan2357]

Fixes #1991

This is a work in progress which implements the bfast algorithm adapted for MEEP's UPML formulation. Will post mathematical notes soon.

To do:

• get it working on a minimal test case
• add tests
• documentation
• remove changes to example files
• convert mathematical documentation to markdown

[codecov-commenter]

Codecov Report

Merging #2609 (b382ffe) into master (cd9a7eb) will increase coverage by 0.16%.
Report is 21 commits behind head on master.
The diff coverage is 100.00%.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2609      +/-   ##
==========================================
+ Coverage   73.90%   74.07%   +0.16%     
==========================================
  Files          18       18              
  Lines        5293     5396     +103     
==========================================
+ Hits         3912     3997      +85     
- Misses       1381     1399      +18     

Files	Coverage Δ
python/simulation.py	77.38% <100.00%> (+0.31%)	⬆️

... and 7 files with indirect coverage changes

[oskooi]

Thanks for working on this feature. One suggestion for a quick test in 1d would be computing the reflectance from a flat interface and verifying the result using the Fresnel equations as demonstrated in this tutorial.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[oskooi]

In case you missed it, there is a type-conversion error in step_db.cpp in the latest commit when running make for this branch during the CI for the single-precision floating-point build as reported in the test logs:

../../../src/step_db.cpp: In member function ‘bool meep::fields_chunk::step_db(meep::field_type)’:
../../../src/step_db.cpp:129:36: error: conversion from ‘vector<double>’ to non-scalar type ‘vector<float>’ requested
  129 |           std::vector<realnum> k = bfast_k_bar;
      |                                    ^~~~~~~~~~~

You can view these test logs by clicking on the Details button for the failing test.

Also, to see the actual error in the failing tests/2D_convergence.cpp, you will need to view the artifacts file i.e. cpp-tests-mpi-true-log by clicking on Summary at the top left of the CI page:

============================================
   meep 1.28.0-beta: tests/test-suite.log
============================================

# TOTAL: 20
# PASS:  19
# SKIP:  0
# XFAIL: 0
# FAIL:  1
# XPASS: 0
# ERROR: 0

.. contents:: :depth: 2

FAIL: 2D_convergence
====================

Using MPI version 3.1, 2 processes
Running holes square-lattice resolution convergence test.
Checking convergence for ey field...
(The correct frequency should be 0.179944.)
frequency for a=30 is 0.179632, 0.179749 (shifted), 0.17969 (mean)
Unshifted freq error is -0.281165/30/30
Shifted freq error is -0.175337/30/30
Frequency difference with a of 30 is -0.105828/30/30
frequency for a=25 is 0.17963, 0.179856 (shifted), 0.179743 (mean)
Unshifted freq error is -0.196313/25/25
Shifted freq error is -0.0551143/25/25
Frequency difference with a of 25 is -0.141199/25/25
frequency for a=20 is 0.179822, 0.180153 (shifted), 0.179988 (mean)
Unshifted freq error is -0.048724/20/20
Shifted freq error is 0.0837453/20/20
Frequency difference with a of 20 is -0.132469/20/20
frequency for a=15 is 0.180183, 0.180115 (shifted), 0.180149 (mean)
Unshifted freq error is 0.0538094/15/15
Shifted freq error is 0.0385003/15/15
Frequency difference with a of 15 is 0.0153091/15/15
frequency for a=10 is 0.180252, 0.181218 (shifted), 0.180735 (mean)
Unshifted freq error is 0.0307579/10/10
Shifted freq error is 0.127421/10/10
Frequency difference with a of 10 is -0.0966632/10/10
frequency for a=5 is 0.181509, 0.187457 (shifted), 0.184483 (mean)
Unshifted freq error is 0.0391222/5/5
Shifted freq error is 0.187825/5/5
Frequency difference with a of 5 is -0.148703/5/5
Passed 2D resolution convergence test for ey!
Checking convergence for ez field...
... using exp(i beta z) z-dependence with beta=0.1
(The correct frequency should be 0.173605.)
meep: Frequency doesn't converge properly with a.
frequency for a=30 is 0.166967, 0.166955 (shifted), 0.166961 (mean)
Unshifted freq error is -5.97385/30/30

[stevengj]

@Dan2357, great to see that you are making progress on this. Let me know if you want to meet sometime.

[stevengj]

@oskooi said he will write the user documentation / tutorial and a test (based on your example — just comparing the reflection to the analytical Fresnel coefficients).

@Dan2357 is going to tweak the technical writeup a bit — expand the intro, and add a section about how this was implemented in Meep as a separate step_bfast function that just adds the new terms, so that we don't need to modify the main timestepping code. He'll also delete the obsolete TeX and PDF files from this PR since everything is now in the .md file.

We should also finalize the API — probably the needs_bfast option can be removed. kbar=None will the be default, and if you pass any other kbar that will enable the bfast option. Not sure if we should call it kbar or something else? (It's just k/ω)

(The technical writeup will then be linked into the main manual somewhere.)

[Dan2357]

I have edited the technical writeup and deleted the unnecessary files.

[oskooi]

edit (11/20): disregard, the fields are not blowing up. The Courant parameter simply needed to be reduced.

As reported in #2722, I think there could be a bug in the current implementation related to not taking into account the index of the source medium when specifying the Bloch-periodic wavevector internally. Note that in the notebook demonstration in this branch, the source medium is air which produces the correct reflectance matching the Fresnel equations whereas for the unit test in #2722 the source medium is $n = 1.4$ and the results do not agree.

As additional information, in the unit test in #2722 I am using:

bfast_k_bar = (np.sin(theta_rad), 0, 0)

which produces the wrong results. If I modify this to include the source medium:

bfast_k_bar = (self.n1 * np.sin(theta_rad), 0, 0)

the fields blow up.

The bug fix therefore likely needs to be made internally rather than by the user when specifying bfast_k_bar

[stevengj]

Note that, in the longer run, it would be good to have a option that lets you put the exp(ikx) factor back into the Fourier-transformed fields, e.g. for mode extraction or visualization. (This can only be done in the DFT fields, not the time-domain fields, because k = k̄ω is ω-dependent.)

Not something for this PR, just something to keep in mind for the future.

[stevengj]

Could we correct the default Courant factor automatically when kbar is specified?

[stevengj]

Let's merge this for now, and maybe update the API in future PR to (1) eliminate needs_bfast, (2) change bfast_kbar to bfast_scaled_k, and (3) automatically update the Courant factor.