[Users] Implementing a simple 1+1 BSSN in Python

Tue May 17 04:01:46 CDT 2022

Hi Ian, Erik,

Thanks for your feedback. I've finally tracked the problem to what I thought the ET was doing at the boundary vs what was actually happening. My par file has

Coordinates::outer_boundary_size = 1

and together with the none BC, the initial values of the system variables remain as such at these points. This is different to my Python code, which is evolving the outer boundary points (in preparation for the SAT method). Thus, I set the time derivative of all fields to zero at the outer boundary points and switched back to RK4 (Euler is indeed unstable!) and everything works now. So, the issue can be summarised as my misunderstanding of what is happening in the ET.

Thanks again for your helpful comments!

Chris

[cid:bea4e148-8810-4b51-bc90-c361db765aeb]

[cid:6b3e7f6f-bae6-4079-9c3d-40fb08baad27]

Dr Chris Stevens

Lecturer in Applied Mathematics

Rm 602, Jack Erskine building

School of Mathematics and Statistics

T: +64 3 369 0396 (Internal 90396)

University of Canterbury | Te Whare Wānanga o Waitaha

Private Bag 4800, Christchurch 8140, New Zealand

http://www.chrisdoesmaths.com<http://www.chrisdoesmaths.com/>

Director
SCRI Ltd
http://www.scri.co.nz<http://www.scri.co.nz/>

________________________________
From: Ian Hinder <ian.hinder at manchester.ac.uk>
Sent: 17 May 2022 19:32
To: Chris Stevens <chris.stevens at canterbury.ac.nz>
Cc: Erik Schnetter <schnetter at gmail.com>; Einstein Toolkit Users <users at einsteintoolkit.org>
Subject: Re: [Users] Implementing a simple 1+1 BSSN in Python

You don't often get email from ian.hinder at manchester.ac.uk. Learn why this is important<https://aka.ms/LearnAboutSenderIdentification>
Hi Chris,

- Note that the Euler method you are using is unstable, and will probably amplify any roundoff errors you have.
- Does everything work if you use periodic boundary conditions?

Erik said:

Given your description, I assume that all derivatives would be identically zero, that there would not even be a floating-point round-off error. Is that so? If not, why not?

I think this would be true for something like (u_{j+1} - u_j) / dx, but it might not be true for ( u_{j+2} - 2u_{j+1} + u_j ) / dx^2 which would evaluate to ( u - 2u + u ) / dx^2, and depending on how the compiler implements it, (u - 2 u) for example might not be exactly representable in floating point, and might not cancel exactly the (+u). If u were zero, I would agree.

On 14 May 2022, at 04:22, Chris Stevens <chris.stevens at canterbury.ac.nz<mailto:chris.stevens at canterbury.ac.nz>> wrote:

Hi Erik,

Thanks for your quick response!

I use the none bc type for now, rather than the more complicated radiative bc. I use a stencil that leans over as you approach the boundary. Currently this is a 4th order interior 3rd order sbp fd operator of Strand.

By looking at the output on the boundary in both codes, there are round off level numbers appearing after the first Euler step in multiple variables. There is no use comparing these exactly as they're so small but note that the same variables have these in both codes. These come from taking spatial derivatives.

Sorry, the lapse is initially one.

I have compared output by saving the data in both codes. It seems like, since the first few time steps everything is round off errors, it is hard to see exactly where things start to go wrong.

I will try your small grid and see if that yields anything more obvious.

Thanks,

Chris

________________________________
From: Erik Schnetter <schnetter at gmail.com<mailto:schnetter at gmail.com>>
Sent: Saturday, May 14, 2022 3:09:17 PM
To: Chris Stevens <chris.stevens at canterbury.ac.nz<mailto:chris.stevens at canterbury.ac.nz>>
Cc: Einstein Toolkit Users <users at einsteintoolkit.org<mailto:users at einsteintoolkit.org>>
Subject: Re: [Users] Implementing a simple 1+1 BSSN in Python

Chris

What do you mean by "Boundary conditions are not applied"? What do you do at the boundaries? Are you using one-sided differencing there?

Given your description, I assume that all derivatives would be identically zero, that there would not even be a floating-point round-off error. Is that so? If not, why not?

Are you really setting the lapse to zero initially? It should be one for Minkowski, same as e.g. g_xx.

What happens if you choose a very small grid size (e.g. just 10 points) and compare the initial conditions between the Einstein Toolkit and the Python code? Are they identical? You would need to compare ASCII output and compare all the digits.

What happens if you take a single Euler time step? Are the values still all identical? Try to find out and which iteration the first grid point differs between the two codes.

-erik

On Fri, May 13, 2022 at 9:59 PM Chris Stevens <chris.stevens at canterbury.ac.nz<mailto:chris.stevens at canterbury.ac.nz>> wrote:
Hi all,

I've been having trouble implementing the BSSN equations, as given by CTGamma, in Python using the package COFFEE, which I am very familiar with:

https://gitlab.com/thebarista/coffee
https://doi.org/10.1016/j.softx.2019.100283

The main comments are:

-- For simplicity, I am starting with a 1+1 code (set 2/3 spatial derivatives to be exactly zero).
-- The 1+log slicing for the lapse and a zero shift.
-- The "phi" conformal factor type
-- The evolution equations are written to Python by a very slight modification of the .m file in CTGEvolution. I've confirmed this .m file outputs the same file that is in the src directory.
-- The vanishing of the trace of A is enforced after each full timestep in the same way as CTGamma.
-- Boundary conditions are not applied.
-- Fourth order SBP operators that lean over toward the boundaries are used.
-- Minkowski initial data with all fields zero except \gamma11/22/33.
-- Euler step is used for time integration to easier compare between Python and the Einstein toolkit.

This leads to a stable evolution in the Einstein toolkit with the attached .par file using a simple one patch Cartesian grid. However, in Python, I get that fields diverge, starting at the boundaries and then propagating inward (including with RK4). This behaviour does not go away with higher resolution or reduced CFL. It is however, fixed completely by setting the second spatial derivative of \alpha to be exactly zero. Trying many different SBP FD operators that have proved successful for other projects has not changed the behaviour.

The COFFEE code has been used and tested in a variety of projects and so there should not be a problem here. Thus, the only thing I can think of, is that the Einstein toolkit is doing something that I am not realizing, that is helping it remain stable. Hopefully, somebody on this mailing list might have an idea?

Thanks in advance,

Chris

<Outlook-zot52qe4.png>

<Outlook-sdnm0er0.png>

Dr Chris Stevens
Lecturer in Applied Mathematics
Rm 602, Jack Erskine building
School of Mathematics and Statistics
T: +64 3 369 0396 (Internal 90396)
University of Canterbury | Te Whare Wānanga o Waitaha
Private Bag 4800, Christchurch 8140, New Zealand
http://www.chrisdoesmaths.com<http://www.chrisdoesmaths.com/>

Director
SCRI Ltd
http://www.scri.co.nz<http://www.scri.co.nz/>

_______________________________________________
Users mailing list
Users at einsteintoolkit.org<mailto:Users at einsteintoolkit.org>
http://lists.einsteintoolkit.org/mailman/listinfo/users

--
Erik Schnetter <schnetter at gmail.com<mailto:schnetter at gmail.com>>
http://www.perimeterinstitute.ca/personal/eschnetter/

_______________________________________________
Users mailing list
Users at einsteintoolkit.org<mailto:Users at einsteintoolkit.org>
http://lists.einsteintoolkit.org/mailman/listinfo/users

--
Ian Hinder
Research Software Engineer
University of Manchester, UK

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.einsteintoolkit.org/pipermail/users/attachments/20220517/abb3176e/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-ldmmb54u.png
Type: image/png
Size: 13436 bytes
Desc: Outlook-ldmmb54u.png
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20220517/abb3176e/attachment-0002.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Outlook-kg41ljz3.png
Type: image/png
Size: 17337 bytes
Desc: Outlook-kg41ljz3.png
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20220517/abb3176e/attachment-0003.png 