Numerical details how to evolve CCZ4 with FD
Collection of what was written by Dumbser in various E-Mails in order to proceed in the Runge Kutta - Finte Differencing code (Cactus/Antelope/Okapi):
Dumbser, 3. April 2018 um 11:03: FO-CCZ4 with Finite differencing works
since Sven has reported some difficulties with the implementation of FO-CCZ4 in the Einstein toolkit last week, and since I wanted to understand the potential problems in depth, I have simply written my own finite difference code for the Einstein equations, based on central finite differences in space and Runge-Kutta time integration. According to Sven's and Elias' description, this is exactly what you are also doing in Cactus, right? The implementation is straightforward, since FD schemes are extremely simple.
To do my tests, I have just copy-pasted my Fortran subroutine PDEFusedSrcNCP into the finite difference code, and I then insert FD point values of Q and central FD approximations for the first spatial derivatives. To save CPU time, I have done all computations in 1D so far.
Please find attached the results that I have obtained for the Gauge wave with amplitude A=0.01 and A=0.1 until a final time of t=1000. I have used sixth order central FD in space and a classical third order Runge-Kutta scheme in time. For FO-CCZ4, I have set all damping coefficients to zero (kappa1=kappa2=kappa3=0), and I use c=0 with e=2. Zero shift (sk=0) and harmonic lapse. CFL number based on e is set to CFL=0.5 at the moment. Now the important points:
Runge-Kutta O3 is for now preferrable over Runge-Kutta O4, since it is intrinsically dissipative. The reason is that the fourth order time derivative term in the Taylor series on the left hand side remains with RK3, while it cancels with RK4, and when moving the term q_tttt to the right hand side and after Cauchy-Kowalevsky procedure, it becomes a fourth order spatial derivative term with negative sign, which is good for stability (second spatial derivatives must have positive sign on the right hand side, fourth spatial derivatives must have negative sign for stability. this is easy to check via Fourier series and the dispersion relation of the PDE).
The RK3 alone is not enough to stabilize the scheme for the larger amplitude A=0.1 of the Gauge Wave, but it is sufficient for A=0.01. I therefore explicitly needed to subtract a term of the type - dt/dx*u_xxxx, which is essentially a fourth order Kreiss-Oliger-type dissipation with appropriately chosen viscosity coefficient.
I will now replace the Kreiss-Oliger dissipation which I do not like with the numerical dissipation that you would have obtained with a Rusanov flux in a fourth order accurate finite volume scheme. In the end, the dissipation operator will again be written as a finite difference, but there I know at least exactly what is going on and we will exactly know the amount of dissipation to be put (it can only be a function of the largest eigenvalue). So there will be NO parameter to be tuned. I will keep you updated on this.
From my own results I can conclude that everything is working as expected, i.e. you must have at least one bug in your implementation of FO-CCZ4 in Cactus; or you have run the tests with the wrong parameters (please use CFL=0.5 for the moment, Kreiss-Oliger dissipation with a viscosity factor so that you get -dt/dx*u_xxxx on the right hand side in the end, please set kappa1=kappa2=kappa3=0 and set e=2 and c=0 in FO-CCZ4). If your code still does not run, I can send you my finite difference Fortran code to help you with the debugging.
While running and implementing my FD code for FO-CCZ4 I have also been working on the vectorization of FO-CCZ4 in ADER-DG. The interesting news: the finite difference code requires more than 20 microseconds per FD point update, and ADER-DG with the good new initial guess for the space-time polynomial and proper vectorization needs also about 20 microseconds per DOF update for FO-CCZ4, i.e. ADER-DG is indeed becoming competitive with FD, who would have every believed this last year :-) On the new 512bit vector machines of RSC in Moscow, we expect the PDE to run even twice as fast, since the vector registers are twice as large as the current state of the art. We are aiming at a time per DOF update of about 10 microseconds. I will keep you informed.
Dumbser, 4. April 2018 um 10:51:
However, my latest experiments show that you can also use RK4 together with a finite-volume type dissipative operator, which is very simple to implement and which does not require any parameters to be tuned. It will just replace the Kreiss-Oliger dissipation. And by the way: in this setting, the scheme can be run with CFL=0.9, which is what we want. I will send around more details later.
Dumbser, 4. April 2018 um 18:23
there are again good news from the finite difference for FO-CCZ4 front. Instead of your classical Kreiss-Oliger dissipation, I suggest to use the following dissipation operator, which should simply be "added" to the time derivatives of
all quantities on the right hand side, i.e.:
dudt(:,i,j,k) = dudt(:,i,j,k) - 1.0/dx(1)* 3.0/256.0* smax * ( -15.0*u(:,i+1,j,k)-u(:,i+3,j,k)-15.0*u(:,i-1,j,k)+6.0*u(:,i+2,j,k)+20.0*u(:,i,j,k)+6.0*u(:,i-2,j,k)-u(:,i-3,j,k) )where dudt(:,i,j,k) is the time derivative of the discrete solution computed by the existing Fortran function PDEFusedSrcNCP and smax is the maximum eigenvalue in absolute value. This operator derives from
- 1/dx(1)*( fp - fm ),where the dissipative flux fp is defined as
fp = - 1/2 * smax * ( uR - uL ),and uR and uL are the central high order polynomial reconstructions of u evaluated at the cell interface x_i+1/2. The flux fm is the same, but on the left interface x_i-1/2.
Please find attached the new results for the Gauge Wave with A=0.1 amplitude. Everything looks fine, i.e., the ADM constraints as well as the waveform at the final time. Note that this simulation was now run with the fourth order Runge-Kutta scheme in time and using a CFL number of CFL=0.9 based on the maximum eigenvalue.