Chris That's good to hear. (I assumed, same as you, that the Einstein Toolkit would evolve the boundary points as well in this case...)
-erik On Tue, May 17, 2022 at 5:01 AM Chris Stevens < chris.stev...@canterbury.ac.nz> wrote: > 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 > > > > > > > > *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 > > *Director* > SCRI Ltd > http://www.scri.co.nz > > ------------------------------ > *From:* Ian Hinder <ian.hin...@manchester.ac.uk> > *Sent:* 17 May 2022 19:32 > *To:* Chris Stevens <chris.stev...@canterbury.ac.nz> > *Cc:* Erik Schnetter <schnet...@gmail.com>; Einstein Toolkit Users < > users@einsteintoolkit.org> > *Subject:* Re: [Users] Implementing a simple 1+1 BSSN in Python > > You don't often get email from ian.hin...@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.stev...@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 <schnet...@gmail.com> > *Sent:* Saturday, May 14, 2022 3:09:17 PM > *To:* Chris Stevens <chris.stev...@canterbury.ac.nz> > *Cc:* Einstein Toolkit Users <users@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.stev...@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 > > *Director* > SCRI Ltd > http://www.scri.co.nz > > _______________________________________________ > Users mailing list > Users@einsteintoolkit.org > http://lists.einsteintoolkit.org/mailman/listinfo/users > > > > -- > Erik Schnetter <schnet...@gmail.com> > http://www.perimeterinstitute.ca/personal/eschnetter/ > > _______________________________________________ > Users mailing list > Users@einsteintoolkit.org > http://lists.einsteintoolkit.org/mailman/listinfo/users > > > -- > Ian Hinder > Research Software Engineer > University of Manchester, UK > > -- Erik Schnetter <schnet...@gmail.com> http://www.perimeterinstitute.ca/personal/eschnetter/
_______________________________________________ Users mailing list Users@einsteintoolkit.org http://lists.einsteintoolkit.org/mailman/listinfo/users
