GRTLCollaboration / GRTeclyn

Port of GRChombo to AMReX - under development!
BSD 3-Clause "New" or "Revised" License
4 stars 2 forks source link

Issue with fourth order cell-centered interpolation #40

Closed mirenradia closed 2 months ago

mirenradia commented 9 months ago

When this code was ported from Chombo to AMReX, Weiqun implemented the following interpolation methods in AMReX in order to improve the accuracy:

We believe GRChombo uses the following instead:

Unfortunately, these seem to have caused a accuracy regression from GRChombo. The GIFs below show a visualization of the CCZ4 variable $\Theta$ (which is a proxy for the constraint violations in the data and should be damped by the constraint damping) on a slice through the equatorial plane of the first few timesteps of a BBH simulation (params_grteclyn.txt for GRTeclyn and params_grchombo.txt for GRChombo). Note that the zoom nor the number of visualized timesteps is not consistent between the GIFs.

Switching back to the methods analogous to the ones in GRChombo (commit a8b0e11) makes GRTeclyn give similar results to GRChombo:

It doesn't seem to be sufficient to just change one of these types of interpolation. In particular, just switching the fine to coarse interpolation to simple averaging (leaving the coarse to fine interpolation as the new fourth order cell-centred interpolation method) improves things but is still significantly worse than GRChombo:

julianakwan commented 5 months ago

I created a separate test case to try to reproduce the problem that solves the Sine Gordon equation: $$\frac{\partial^2 \phi}{\partial t^2} = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} - \sin(\phi)$$ This exists as separate AMReX standalone code and also within the GRTeclyn framework in the Klein-Gordon branch. The idea here is that we have a lightweight code that runs quickly with an analytic solution that we can check against.

1D Sine Gordon

In the 1D case, an analytic breather solution exists: $$\phi (x,t) = 4 \arctan \left(\frac{\beta}{\alpha}\frac{\cos(\alpha t)}{\cosh (\beta x)}\right),\; \alpha^2 + \beta^2 = 1 $$ which is supposed to continue oscillating in place indefinitely while changing in sign. In the following figures, the error estimates are made with respect to the above solution. In this test, I use the solution as the initial condition with $\alpha = 0.7$ and $t_0= 0$ and allowed for 3 levels of refinement (so max_level = 3). Switching between quartic and cell quartic (by changing the interpolator in when defining desc_lst/deriv_lst and fill_patch/average_down) showed that the cell quartic interpolator was actually more accurate.

Here is a movie of the difference between $\phi$ calculated by AMReX using the quartic interpolator and the exact solution: diff_3_levels_quartic_upload.webm

Same error calculation but between $\phi$ using the cell quartic interpolator and the exact solution: diff_3_levels_cell_quartic_upload.webm

You can see that the maximum difference for the quartic case is ~5 x 10^-3 while for the cell quartic case it is ~10^-3.

3D Sine Gordon

In the 3D case, there is no exact solution and only a pseudo breather: $$\phi(x,y,z,t) = 64 \arctan \left(\frac{\beta}{\alpha}\frac{\cos(\alpha t_0)}{\cosh (\beta x)}\right)\arctan\left(\frac{\beta}{\alpha}\frac{\cos(\alpha t_0)}{\cosh (\beta y)}\right)\arctan\left(\frac{\beta}{\alpha}\frac{\cos(\alpha t_0)}{\cosh (\beta z)}\right)$$ I used $\alpha = 0.7$ and two blobs: one at (50,50,25) with $t_0 = -5.7$ and the other at (50,50,75) with $t_0 = +5.7$. The full extent of the grid is (100,100,100). (The second blob was for testing the ray-tracing algorithm, Wombat)

Because there is no exact solution, the error estimates are made with respect to a run with no AMR and the entire region is covered by a grid with the same spacing as the finest level. In this case, the two test runs (using the quartic and cell quartic interpolators) have max_level = 2 with N_full = 64 so the finest level is N = 256. The "coarse grid" only run for comparison purposes is then N_full = 256. I used an initial dt = 0.25 for all three. I then excised the level 2 grids from each of the cell quartic and quartic runs and compared the values of $\phi$ to the corresponding region in the no AMR run. Again, I found that the run with the cell quartic interpolator were more accurate: it had a lower mean squared error measured per cell (0.0057 vs 0.0102) and a lower median fractional error (0.0104 vs. 0.1777). I made a histogram of the fractional error per cell below: sine_gordon_3d_frac_error_histogram.pdf You can see that, in both cases, 48 cells have some really extreme values (from the ghost regions?), but the fractional error for the cell quartic interpolator (orange) is always lower than that for the quartic interpolator (blue).

mirenradia commented 5 months ago

Following a discussion with @KAClough, I experimented with replacing the explicit trace removals of $\tilde{A}_{ij}$ with the addition of conformal constraint damping terms as described in Ref. ^1. This can be found on the test/conformal_constraint_damping branch. Unfortunately, this didn't seem to help as can be seen from the visualization of $\Theta$ (commit f9b0b989878bb5c17bed564311d07c351fadd23a): GRTeclyn-Theta-conformal-constraint-damping

KAClough commented 5 months ago

(There is no crying emoji)

KAClough commented 4 months ago

I noted an error in the implementation of the dynamic trace removal. Having fixed it in this branch https://github.com/GRTLCollaboration/GRTeclyn/tree/test/conformal_constraint_damping and tested some different values of kappac I can confirm that it still makes no difference to the behaviour 😢 (but in other better news, I found the crying emoji)

mirenradia commented 4 months ago

I noted an error in the implementation of the dynamic trace removal. Having fixed it in this branch https://github.com/GRTLCollaboration/GRTeclyn/tree/test/conformal_constraint_damping and tested some different values of kappac I can confirm that it still makes no difference to the behaviour 😢 (but in other better news, I found the crying emoji)

@KAClough, for completeness, would you be able to push a commit with the fixed code to that branch?

KAClough commented 4 months ago

Yes, I think I already did this here https://github.com/GRTLCollaboration/GRTeclyn/commit/b735f2e4528538f840215dd2eb9862c57b19107e

mirenradia commented 4 months ago

Ah, my bad. For some reason, I was expecting a notification on the gitfeed on Slack but these aren't enabled for non-default branches at the moment.

julianakwan commented 3 months ago

I made a quick video of the BinaryBH example based on the bugfix40 branch:

https://github.com/GRTLCollaboration/GRTeclyn/assets/27848011/f9f8bb25-4757-4b87-8cdf-fcc230c223cb

(NB: this is only 85 timesteps)

mirenradia commented 2 months ago

For completeness, here is the same animation of $\Theta$ in the the first comment after the fix in 0a892700a6941cb2eae6ca80c92c1a1541cde4c6:

GRTeclyn-Theta-fixed

I am confident this issue is resolved.