afd-illinois/grim

General Relativistic Implicit Magnetohydrodynamics

grim's Issues

Possible errors in fluidElement due to redundant input

fluidElement constructor takes const grid &prim and const geometry &geom as inputs, but its functions also take these objects as inputs. This can lead to inconsistencies and errors if, for example, a fluidElement is constructed with centered geometry and primitive variables but its functions are called with geometry and primitive variables evaluated on cell faces.

primOldLocal not set during second half step in EXPLICIT and IMEX evolutions

Bug discovered by @ffoucart

During the second half step, we have the flags:
computeDivFluxAtTimeN = 0
computeDivFluxAtTimeNPlusHalf = 1

The above configuration sets primHalfStepLocal in the residual computation and primOldLocal is not set. However, primOldLocal is still being used to computed conservedOldVars.

Fix: Need to change primOldLocal to primOldGlobal before computing conservedOldVars.

Also add a flag to avoid computing conservedOldVars during the second half step since they are already computed during the first half step

Issues with linear modes in the conduction code

1. Convergence is at 1.7 order instead of 2nd order. Probably because I'm not comparing the numerical solution and the theoretical solution at the exact same time. O(dt) error? The spatial derivatives seem to converge at 2nd order as seen from the atmosphere test.

2. Why are local extrema being produced? A possible bug or is there any better way to compute derivatives in computeConductionSpatialGradientTerms() ?
   

Configure fails with latest version of petsc-dev

Use commit fced3c3f9e703542693913793d15321603e40fe6 while compiling petsc-dev in the meantime.

Should functions take grids or arrays as arguments

Uncertain how to decide between grid and array types for function arguments. Consider for example a function that transforms a four-vector in the coordinate frame vCoord to the tetrad frame vTetrad; should this function be

void fluidElement::coordToTetrad(const array vCoord[NDIM], array vTetrad[NDIM])

or

void fluidElement::coordToTetrad(const grid &vCoord, grid &vTetrad)

?

Wallclock termination requires changes to params files of all problems

The newly implemented method to terminate the code after a given time, and allow checkpoint/restart without recompiling, will break existing problems unless new parameters are added to the params.input files:

//Name of files to store last checkpoint information
std::string restartFileName = "restartFileName.txt";
std::string restartFileTime = "restartFileTime.txt";

as well as the function

int timestepper::CheckWallClockTermination()

which can either return 0 (never terminate due to wallclock), or be copied from the torus.

Anisotropic Israel-Stewart viscosity

Guidelines to start off with:

1. Start with a new branch from master
   git checkout -b viscosity

2. Make data structures in timestepper.h that will hold the spatial gradients of the terms in the viscosity theory. See the code within #IF (CONDUCTION) #ENDIF for reference. Write the init and destroy routines for these data structures.

3. Add the new degree of freedom for the scalar shi in physics.h. Here one has to be careful. We want to be able to run with/without conduction and viscosity independently and so the macro DOF needs to be defined accordingly.

4. Copy physics/conduction.c to physics/viscosity.c. Like conduction.c, viscosity.c will consist of two functions:
   void addViscositySourceTermsToResidual(...)
   void computeViscositySpatialGradientTerms(...)

5. Add the above function in the appropriate locations in timestepper/residual.c. These should be put in the same locations as the conduction calls and with #IF (VISCOSITY) #ENDIF.

6. Add viscosity.c to physics/CMakeLists.txt

7. Add the necessary options such as set(VISCOSITY "ON") to the problem CMakeLists.txt. Also need to add a function in problem/problem.h which will set the transport coefficients for viscosity, just like in conduction. The function definition is to be given problem/whicheverproblem/problem.c

8. Commit changes and push to remote
   git commit -a
   git push -u origin viscosity

Definitions of q and dP

Need to decide what to use as primitive variables and conserved variables.

Strange segfaults on bh27

grim randomly segfaults on bh27 sometimes. Debugger indicates that it has something to do with OpenMP. I thought it is because the tile size is too large but that ain't it.

Bug in time step

computeExplicitSources() being called before computeEMHDGradients() (bug introduced in f5c36cc)

Fixed in 68c9814

Issue opened to notify all clones/forks.

Clean up of the local linesearch code

timeStep() in src/timestepper/timestepper.c is very messy right now.

Preferred approach is to move the local linesearch into a new function and put this function in a new file locallinesearch.c, since there is substantial code in there.

OpenMP doesn't work in 1D

Remember to set USE_OPENMP to NO when using COMPUTE_DIM = 1.

Automated tests

Implement automated tests that can be executed in under a minute that will ensure the correctness of the solutions produced by the code at every commit.

Ideas off the top of my head

SR:

1. Linear modes : The code runs with a known mode for a few wave periods and the computes the norm of the error between the evolved solution and the analytic solution.

GR:
2) Static atmosphere: Evolve a steady state atmosphere which is in pressure equilibrium and compute the norm of the error between the evolved solution and the initial state.

3. Bondi inflow: Same procedure as the static atmosphere case but with the Bondi inflow solution. Can check the values of heat flux and pressure anistropy by running without backreaction.

Suggestions welcome.

Make the code insensitive to checkpoint/restart

When restarting from a dumped file, the code doesn't know what the current time step is. So the restarted evolution is not identical to the original evolution. This can make debugging tricky.

COMPUTE_DIM=1 doesn't work with KERRSCHILD

Atmosphere and Bondi tests fail when used with COMPUTE_DIM=1. Just use COMPUTE_DIM=2 for now.

Copy all problem options in run directory

Put in wavespeeds for conduction and viscosity in the Riemann solver

Seg fault with local linesearch

Hi Francois,

I'm trying to run linear modes with the local linesearch method, but the code segfaults at

(gdb) bt
#0  0x000000000049779c in timeStep (ts=ts@entry=0x7fffffffdce0) at /home/astro/work/grim/src/timestepper/timestepper.c:591
#1  0x0000000000496b5a in main (argc=11, argv=0x7fffffffdee8) at /home/astro/work/grim/src/grim.c:16

I'm tried with all three timestepping options EXPLICIT, IMEX and IMPLICIT, and they all give the errors.

I ran with the code with

./grim -snes_max_it 1 -pc_type asm -ksp_type preonly -mat_fd_type ds -snes_linesearch_type basic

Do you see something obviously wrong in that line of code?

Could you encapsulate the linesearch method in timestepper.c into a function, so that it is cleaner and easier to debug?

Thanks,
Mani

Better problem-dependent initial guess for failed points

Mac compilation

For some reason Petsc's SETERRQ() function doesn't seem to compile on a mac. I think the compiler being used is clang.

Residual normalization

Residuals in GRIM need to be normalized. Right now we only divide by a factor a sqrt(-g).

geom tests failing

The following tests FAILED:
     19 - Minkowski_params_1D_4_procs (Failed)
     20 - Minkowski_gCov_1D_4_procs (Failed)
     21 - Minkowski_gCon_1D_4_procs (Failed)
     22 - Minkowski_g_1D_4_procs (Failed)
     23 - Minkowski_alpha_1D_4_procs (Failed)
     24 - Minkowski_params_2D_4_procs (Failed)
     25 - Minkowski_gCov_2D_4_procs (Failed)
     26 - Minkowski_gCon_2D_4_procs (Failed)
     27 - Minkowski_g_2D_4_procs (Failed)
     28 - Minkowski_alpha_2D_4_procs (Failed)
     29 - Minkowski_params_3D_4_procs (Failed)
     30 - Minkowski_gCov_3D_4_procs (Failed)
     31 - Minkowski_gCon_3D_4_procs (Failed)
     32 - Minkowski_g_3D_4_procs (Failed)
     33 - Minkowski_alpha_3D_4_procs (Failed)
     34 - Kerr_params_1D_4_procs (Failed)
     35 - Kerr_gCov_1D_4_procs (Failed)
     36 - Kerr_gCon_1D_4_procs (Failed)
     37 - Kerr_g_1D_4_procs (Failed)
     38 - Kerr_alpha_1D_4_procs (Failed)
     39 - Kerr_connection_coeffs_1D_4_procs (Failed)
     40 - Kerr_params_2D_4_procs (Failed)
     41 - Kerr_gCov_2D_4_procs (Failed)
     42 - Kerr_gCon_2D_4_procs (Failed)
     43 - Kerr_g_2D_4_procs (Failed)
     44 - Kerr_alpha_2D_4_procs (Failed)
     45 - Kerr_connection_coeffs_2D_4_procs (Failed)
     46 - Kerr_params_3D_4_procs (Failed)
     47 - Kerr_gCov_3D_4_procs (Failed)
     48 - Kerr_gCon_3D_4_procs (Failed)
     49 - Kerr_g_3D_4_procs (Failed)
     50 - Kerr_alpha_3D_4_procs (Failed)
     51 - Kerr_connection_coeffs_3D_4_procs (Failed)
Errors while running CTest
make: *** [test] Error 8

Should the primitive variables be floored inside the nonlinear solver?

Note: If so, then don't exceed 1e-15, else SNES will NOT converge.