Add special relativistic solver (#86)

This adds a special relativistic solver. It includes three Riemann solvers (CGF-style solvers based on the conservative and primitive variables, plus an HLLC solver), tests and problem setups.

* Gresho vortex for lm

* compressible gresho vortex

* SR compressible

* sr advection works??

* corrected RHLLC solver so advects in the right direction now

* put c2p routines in a fortran function so it's no so slow

* fortran c2p is much much faster

* advect, kh, quad, logo, acoustic pulse all seem to work for compressible_sr. gresho, sod and anything with gravity do not - may be due to dodgey initial data or more serious

* Fixed tests for compressible_sr

* Added scipy to requirements file

* Made flake8 compliant

* p != h

* Added compressible_sr example to solver-test notebook

* Updated h5py so no longer get error warning in first cell of solver-test

* Added codecov

* Added more stuff to solver-test so will test more problems and more of compressible_sr/derives.py

* Fixed qy/qx in c2p_f so now sod works

* shhh flake8

* rt and rt2 work!

* bubble and sedov work

* All problems now work?#

* Fixed patch.py for external boundaries

* shhh flake8

* shhh flake8

* Got CGF solver working for SR compressible

* eigenvalues in cgf prim

* updating docs

* updating docs

* Switched to reconstruct using the conservative variables as cannot find full primitive SRHD eigenstructure

* flake8

* Fixed bubble(?)

* converted c2p_f.f90

* converted inteface_f.f90 to python

* make unsplit fluxes compatible with python interface

* fixed a couple of bugs so compressible_sr seems to work now

* delete Fortran files

Author: harpolea

.travis.yml

@@ -26,6 +26,12 @@ before_script:
 script:
   - flake8 .
   - pytest -v --cov=. --cov-config .coveragerc --nbval --ignore=docs --ignore=./multigrid/variable_coeff_elliptic.ipynb --ignore=examples/mesh --ignore=examples/multigrid
+
+after_success:
+    - codecov
+after_failure:
+   - codecov
+   
 #   - travis-sphinx build
 
 # after_success:

compressible/problems/_gresho.defaults

@@ -0,0 +1,9 @@
+[gresho]
+dens_base = 10.0            ; density at the base of the atmosphere
+scale_height = 2.0          ; scale height of the isothermal atmosphere
+
+r = 1.0
+u0 = 1.0
+p0 = 1.0
+
+dens_cutoff  = 0.01

compressible/problems/gresho.py

@@ -0,0 +1,81 @@
+from __future__ import print_function
+
+import sys
+import numpy
+import mesh.patch as patch
+from util import msg
+
+
+def init_data(my_data, rp):
+    """ initialize the Gresho vortex problem """
+
+    msg.bold("initializing the Gresho vortex problem...")
+
+    # make sure that we are passed a valid patch object
+    if not isinstance(my_data, patch.CellCenterData2d):
+        print("ErrrrOrr: patch invalid in bubble.py")
+        print(my_data.__class__)
+        sys.exit()
+
+    # get the density and velocities
+    dens = my_data.get_var("density")
+    xmom = my_data.get_var("x-momentum")
+    ymom = my_data.get_var("y-momentum")
+    ener = my_data.get_var("energy")
+
+    # grav = rp.get_param("compressible.grav")
+
+    gamma = rp.get_param("eos.gamma")
+
+    # scale_height = rp.get_param("gresho.scale_height")
+    dens_base = rp.get_param("gresho.dens_base")
+    dens_cutoff = rp.get_param("gresho.dens_cutoff")
+
+    rr = rp.get_param("gresho.r")
+    u0 = rp.get_param("gresho.u0")
+    p0 = rp.get_param("gresho.p0")
+
+    # initialize the components -- we'll get a psure too
+    # but that is used only to initialize the base state
+    xmom[:, :] = 0.0
+    ymom[:, :] = 0.0
+    dens[:, :] = dens_cutoff
+
+    # set the density to be stratified in the y-direction
+    myg = my_data.grid
+    pres = myg.scratch_array()
+
+    dens[:, :] = dens_base
+
+    pres[:, :] = p0
+
+    x_centre = 0.5 * (myg.x[0] + myg.x[-1])
+    y_centre = 0.5 * (myg.y[0] + myg.y[-1])
+
+    rad = numpy.sqrt((myg.x2d - x_centre)**2 + (myg.y2d - y_centre)**2)
+
+    pres[rad <= rr] += 0.5 * (u0 * rad[rad <= rr]/rr)**2
+    pres[(rad > rr) & (rad <= 2*rr)] += \
+        u0**2 * (0.5 * (rad[(rad > rr) & (rad <= 2*rr)] / rr)**2 +
+        4 * (1 - rad[(rad > rr) & (rad <= 2*rr)]/rr +
+        numpy.log(rad[(rad > rr) & (rad <= 2*rr)]/rr)))
+    pres[rad > 2*rr] += u0**2 * (4 * numpy.log(2) - 2)
+    #
+    uphi = numpy.zeros_like(pres)
+    uphi[rad <= rr] = u0 * rad[rad <= rr]/rr
+    uphi[(rad > rr) & (rad <= 2*rr)] = u0 * (2 - rad[(rad > rr) & (rad <= 2*rr)]/rr)
+
+    xmom[:, :] = -dens[:, :] * uphi[:, :] * (myg.y2d - y_centre) / rad[:, :]
+    ymom[:, :] = dens[:, :] * uphi[:, :] * (myg.x2d - x_centre) / rad[:, :]
+
+    ener[:, :] = pres[:, :]/(gamma - 1.0) + \
+                0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]
+
+    eint = pres[:, :]/(gamma - 1.0)
+
+    dens[:, :] = pres[:, :]/(eint[:, :]*(gamma - 1.0))
+
+
+def finalize():
+    """ print out any information to the userad at the end of the run """
+    pass

compressible/problems/inputs.gresho

@@ -0,0 +1,33 @@
+# simple inputs files for the unsplit CTU hydro scheme
+
+[driver]
+max_steps = 2000
+tmax = 10000.0
+cfl = 0.8
+
+
+[io]
+basename = lm_gresho_128_
+n_out = 1
+
+
+[mesh]
+nx = 128
+ny = 128
+xmax = 1.0
+ymax = 1.0
+
+xlboundary = periodic
+xrboundary = periodic
+
+ylboundary = periodic
+yrboundary = periodic
+
+
+[gresho]
+r = 0.2
+u0 = 0.1
+p0 = 1
+
+[compressible]
+grav = 0

compressible/unsplit_fluxes.py

@@ -316,7 +316,7 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt):
 
     """
     finally, we can construct the state perpendicular to the interface
-    by adding the central difference part to the trasverse flux
+    by adding the central difference part to the transverse flux
     difference.
 
     The states that we represent by indices i,j are shown below