Refactor GRHD project (!101)

* update examples once more

* remove comments

* cleanup: remove codegen for now

* IntegartionTests: add grhd_tov_spherical1d ref test

* IntegrationTests: remove reltol from grhd_bondi_accretion test

* spherical1d: fix missing shift from maxspeed computation

* cleanup

* make spherical1d formulation work when r=0 is included by removing potential nans in source terms by enforcing symmetry

* make impose_symmetry! dispatch on project

* add tov example that runs with spherical1d formulation

* wip: change maxspeed implementation for tov star initial data

* cons2prim: use non-zero values for rho, eps thresholds

* fix: set atm density in bondi_accretion initial data

* implement tov initial data for spherical1d formulation

* add formulation parameter

* implement rhs for FVElement and rescaled_spherical1d

* implement equations for FVElement and rescaled_spherical1d

* mesh: implement differentiate for FVElement

* add grhd example par files

* update notebooks

* fix missing lapse factor in source terms

* fix typos

* fixup initialdata setup

* derivative computation needs jacobi

* rename cons2prim routines

* add atmosphering to cons2prim

* implement rhs! for rescaled_spherical1d formulation

* specialize timespeed computation based on equations

* add specialized equations for the rescaled GRHD version

* implement initialdata setup for tov case

* spherical accretion:make sure K, Gamma agree with loaded data

* expose spectral differentiation

* wrap SphericalAccretion solver into module; recompute bondi_accretion.h5 with less radial points

* refactor TOV solver to use OrdinaryDiffEq

* initialdata: cleanup dir

* add bondi accretion as ref test

* fix up loading of id file

* fv_rhs: implement fv_update_step! with source term

* add ROOTDIR const global so that parfiles can specify paths relative to project root

* update SphericalAccretion.jl

* update sympy nb

* cleanup initialdata setup

* equations: fix up Tmunu computation

* Revert "fix source term in dg and fv rhs computation"

This reverts commit 161cc034.

* Replace Interpolations.jl with BasicInterpolators.jl

* equation: add missing flux terms to sources; use outflow bdry condition on inner radius

* comments in cons2prim

* implement momentum constraint for debugging

* test SphericalAccretion

- copy Bugner's C version of the solver
- compare my results with Bugner's -- good agreement
- compare against literature results (Papadopoulos 1998,1999) --
  discrepancy in sign of vr, ur, otherwise good agreement

* SphericalAccretion: fix wrong definition of uc^r and Vc

* tweak impl more

* fix source term in dg and fv rhs computation

* fix equations and implement fv method

* equations: combine flux and source computation

* math: update spherical_accretion notebook

* SRHD: export module

* re-implemented GRHD eqs with a non-zero shift

* initialdata: add grhd_source.codegen.jl that was generated with TensorComponents.jl

* deps: need Interpolations.jl for initial data setup

* update math/spherical_accretion.py

* update setup

* start initialdata setup for bondi problem

* rework rhs.jl

* update deps and parameters

* update type defs

* initialdata: add HDF5 dep; write h5 output for SphericalAccretion.jl

* rename cons2prim.jl

* remove unused files

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.0-rc2"
 manifest_format = "2.0"
-project_hash = "4772e5e6caf338d2c755b5d4461c2d0c5c7482c5"
+project_hash = "7e9cb4087253abe43f5a39d1ae318381c828bb4b"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
@@ -46,6 +46,12 @@ uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[deps.BasicInterpolators]]
+deps = ["LinearAlgebra", "Memoize", "Random"]
+git-tree-sha1 = "3f7be532673fc4a22825e7884e9e0e876236b12a"
+uuid = "26cce99e-4866-4b6d-ab74-862489e035e0"
+version = "0.7.1"
+
 [[deps.BitTwiddlingConvenienceFunctions]]
 deps = ["Static"]
 git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b"
@@ -409,6 +415,12 @@ deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 version = "2.28.2+1"
 
+[[deps.Memoize]]
+deps = ["MacroTools"]
+git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
+uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
+version = "0.4.4"
+
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "b01beb91d20b0d1312a9471a36017b5b339d26de"

Project.toml

@@ -5,6 +5,7 @@ version = "1.2.0"
 
 [deps]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+BasicInterpolators = "26cce99e-4866-4b6d-ab74-862489e035e0"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Jacobi = "83f21c0b-4282-5fbc-9e3f-f6da3d2e584c"
@@ -27,12 +28,12 @@ ToggleableAsserts = "07ecc579-1b30-493c-b961-3180daf6e3ae"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
-julia = "1.9"
 PrettyTables = "<2.3.0"
+julia = "1.9"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","InteractiveUtils"]
+test = ["Test", "InteractiveUtils"]

examples/grhd_bondi/dg.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 0.02143872868864732
+polytrope_gamma = "$(4/3)"
+
+[GRHD]
+bc = "from_id"
+id = "bondi_accretion"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"bondi_accretion.h5\"))"
+# formulation = "spherical1d"
+
+[Mesh]
+range = [ 1.8, 10.0 ]
+n = 4
+k = 35
+basis = "lgl"
+periodic = false
+
+[Output]
+every_iteration = 1
+# every_dt  = 0.1
+variables = [ "D", "Sr", "τ", "p", "ρ", "ϵ", "vr", "rhs_D", "rhs_Sr", "rhs_τ", "src_D", "src_Sr", "src_τ" ]
+# aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d  = true
+
+[Evolution]
+cfl = 0.4
+tend = 20.0

examples/grhd_bondi/fv.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 0.02143872868864732
+polytrope_gamma = "$(4/3)"
+
+[GRHD]
+bc = "from_id"
+id = "bondi_accretion"
+id_filename = """$(joinpath(ROOTDIR,"initialdata","bondi_accretion.h5"))"""
+
+[Mesh]
+range = [ 1.8, 10.0 ]
+k = 1024
+scheme = "FV"
+periodic = false
+
+[Output]
+every_iteration = 1
+variables       = [ "D", "Sr", "τ", "p", "ρ", "ϵ", "vr", "rhs_D", "rhs_Sr", "rhs_τ", "src_D", "src_Sr", "src_τ" ]
+# aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d        = true
+
+[Evolution]
+cfl = 1.0
+tend = 20.0

examples/grhd_tov/dg.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+
+[Mesh]
+range = [ 0.5, 20.0 ]
+n = 1
+k = 2048
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+every_dt  = 0.2
+variables = [ "D", "Sr", "τ",
+              "p", "ρ", "ϵ", "vr",
+              # "rD", "rSr", "rτ",
+              # "rhs_rD", "rhs_rSr", "rhs_rτ",
+              # "src_rD", "src_rSr", "src_rτ",
+              # "flr_rD", "flr_rSr", "flr_rτ",
+              # "A", "∂Ar", "B", "∂Br"
+              ]
+# aligned_ts      = "$(collect(range(0.01,0.6,step=0.01)))"
+enable1d  = true
+
+[Evolution]
+cfl = 0.2
+tend = 40.0

examples/grhd_tov/dg_spherical1d.toml

+[EquationOfState]
+eos = "polytrope"
+polytrope_k = 100.0
+polytrope_gamma = 2.0
+
+[GRHD]
+bc = "from_id"
+id = "tov"
+id_filename = "$(joinpath(ROOTDIR,\"initialdata\",\"TOV_stable.h5\"))"
+atm_factor = 1e-8
+formulation = "spherical1d"
+
+[Mesh]
+range = [ 0.5, 20.0 ]
+# n = 1
+# k = 2048
+n = 4
+k = 124
+basis = "lgl"
+periodic = false
+
+[Output]
+# every_iteration = 1
+# every_dt  = 0.2
+variables = [ "D", "Sr", "τ",
+              "p", "ρ", "ϵ", "vr",
+              "src_D", "src_Sr", "src_τ",
+              "flr_D", "flr_Sr", "flr_τ",
+              # "A", "∂Ar", "B", "∂Br"
+              ]
+aligned_ts      = "$(collect(range(0.2,120.0,step=0.2)))"
+enable1d  = true
+
+[Evolution]
+cfl = 0.2
+tend = 120.0

initialdata/.gitignore

+SphericalAccretion
+*.txt

initialdata/Makefile

+default:
+	gcc SphericalAccretion.c -o SphericalAccretion -lm
+
+clean:
+	rm SphericalAccretion bondi-acc-ref.txt

initialdata/Project.toml

 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-VimBindings = "51b3953f-5e5d-4a6b-bd62-c64b6fa1518a"
-dg1d = "dfc8a8d9-21b0-4914-8ba8-32d2ec29f3fd"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

initialdata/SphericalAccretion.c

+// copied from bamps/src/projects/grhd/initialdata/michel.c
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+void main(int argc, char *argv[]) {
+
+    double M = 1.0;
+    double rc = 200.0*M;
+    double rhoc = 7e-4;
+    double gamma = 4.0/3.0;
+
+    double urc   = sqrt(.5 * M / rc);
+    double Vc2   = urc*urc / (1.-3.*urc*urc);
+
+    double K     = Vc2*(gamma-1.)*pow(rhoc,1.-gamma) / (gamma*(-Vc2+gamma-1.));
+
+    double hc    = 1. + gamma / (gamma-1.) * K * pow(rhoc,gamma-1.);
+    double c1    = rhoc * urc * rc * rc;
+    double c3    = hc * hc * (urc*urc + 1. - 2.*M / rc);
+
+    int nnn = 1000;
+    int i;
+    double r;
+    double rmin = 1.8*M, rmax = 200.0*M;
+    double dr = (rmax-rmin)/(nnn-1);
+
+    double *gr      = malloc(nnn*sizeof(double));
+    double *grho    = malloc(nnn*sizeof(double));
+    double *geps    = malloc(nnn*sizeof(double));
+    double *gp      = malloc(nnn*sizeof(double));
+    double *gvr     = malloc(nnn*sizeof(double));
+    double *gur     = malloc(nnn*sizeof(double));
+    double *gut2    = malloc(nnn*sizeof(double));
+
+    for(i=0; i<nnn; i++) {
+
+      r = rmax - dr*i;
+      gr[i] = r;
+
+      double ur    = sqrt(fabs(-1. + 2.*M / r) + 0.01);
+      double ut2, h, f, df;
+
+      int it=0;
+      for(it=0;it<100;it++) {
+
+        ut2 = (ur*ur+1.-2.*M/r);
+        h   = 1. + (gamma / (gamma-1.)) * K * pow(c1 / (ur*r*r), gamma-1.);
+        f   = h*h*ut2 - c3;
+        df = 2. * h * (-(h-1.)*(gamma-1.)/ur) * ut2 + 2. * h * h * ur;
+
+        ur = ur - f/df;
+
+        if(fabs(f/df)<1.e-14) break;
+
+      }
+
+      grho[i] = c1 / (ur*r*r);
+      gp[i]   = K * pow(grho[i],gamma);
+      geps[i] = gp[i] / (grho[i]*(gamma-1.));
+
+      ut2 = (ur*ur+1.-2.*M/r);
+      gur[i] = ur;
+      gut2[i] = ut2;
+
+      double alpha = sqrt(r / (r + 2.*M));
+      double betar = 2.*M / (r + 2.*M);
+      double ut    = (-2.*M * ur / r + sqrt(ut2)) / (-1.+2.*M/r);
+      double vr    = ur / (alpha*ut) + betar / alpha;
+
+      gvr[i] = vr;
+
+    }
+
+
+    FILE *io = fopen("bondi-acc-ref.txt", "w");
+    if (!io) {
+      printf("FAILED\n");
+      exit(-1);
+    }
+    for (i=nnn-1; i>=0; i--) {
+      fprintf(io, "% .6e % .6e % .6e % .6e % .6e % .6e % .6e\n", gr[i], grho[i], geps[i], gp[i], gvr[i], gur[i], gut2[i]);
+    }
+    fclose(io); io = NULL;
+
+    free(grho); grho = NULL;
+    free(geps); geps = NULL;
+    free(gp);   gp = NULL;
+    free(gvr);  gvr = NULL;
+    free(gr);   gr = NULL;
+    free(gur);  gur = NULL;
+    free(gut2); gut2 = NULL;
+
+}

initialdata/SphericalAccretion.jl

+module SphericalAccretion
+
+
+
 using CairoMakie
+using CairoMakie.FileIO
+using DelimitedFiles
 using Printf
+using HDF5
+using Roots
+using UnPack
 
 
 """
@@ -7,53 +16,51 @@ Solution taken from Bugner [2017], section 4.2.
 
 Kerr-Schild coordinates allow use to go beyond horizon towards the singularity.
 However, it becomes more difficult for the Newton Rapshon solver to determine a solution
-once beyond the horizon. To this end, try increasing `N_rs` whenever you lower `rmin` and 
+once beyond the horizon. To this end, try increasing `N_rs` whenever you lower `rmin` and
 the solver aborts inside the horizon.
 Alternatively, one could also add an adaptive stepper.
 
 Note: We implicitly assume θ = π/2, because then sin(θ) = 1, although,
 the sines should cancel out everywhere.
 """
-function solve()
-
-  # parameters
-  M = 1.0
-  N_rs = 100000
-  rmin, rmax = 0.1*M, 200*M
-  rc, ρc, Γ = 200*M, 7e-4, 4/3
+function solve(; rmin, rmax, rc, ρc, Γ, M=1.0, N_rs=1000,
+                 urc_sign=1)
 
-  urc = sqrt(M / (2 * rc))
+  urc = urc_sign*sqrt(M /(2 * rc))
   # cannot find solutions within these regions, cf. Michel 1972
   @show urc
   @assert rc >= 6*M
-  @assert urc <= 1 / 3
+  @assert urc^2 <= 1 / 3
 
-  Vc2 = urc^2 / (1 - 3 * urc^2)
-  K = Vc2 * (Γ - 1) / ( (Γ - 1 - Vc2) * Γ * ρc^(Γ - 1) )
+  Vc = urc^2 / (1 - 3 * urc^2)
+  K = Vc * (Γ - 1) / ( (Γ - 1 - Vc) * Γ * ρc^(Γ - 1) )
+  @show K
 
   # here we set θ = π/2 so that sin(θ) = 1
-  C1 = rc^2 * ρc * urc
   hc = 1 + K * Γ * ρc^(Γ - 1) / (Γ - 1)
   u0c = sqrt(1 - 2 * M / rc + urc^2)
+  C1 = rc^2 * ρc * urc
+  C2 = rc^2 * ρc * hc * u0c * urc
   C3 = hc^2 * u0c^2
 
   function F(ur, r)
     ρ = C1 / (r^2 * ur)
     h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
     u0 = sqrt(1 - 2 * M / r + ur^2)
-    u0 = sqrt(1 - 2 * M / r + ur^2)
     return h^2 * u0^2 - C3
   end
 
-  function dF(ur, r)
-    ρ = C1 / (r^2 * ur)
-    h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
-    u0 = sqrt(1 - 2 * M / r + ur^2)
-    dρdur = - C1 / (r^2 * ur^2)
-    dhdρ = K * Γ * ρ^(Γ - 2)
-    du0dur = ur / u0
-    return 2 * (u0^2 * h * dhdρ * dρdur + h^2 * u0 * du0dur)
-  end
+  # manually compute derviative and
+  # use Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton()) below
+  # function dF(ur, r)
+  #   ρ = C1 / (r^2 * ur)
+  #   h = 1 + K * Γ * ρ^(Γ - 1) / (Γ - 1)
+  #   u0 = sqrt(1 - 2 * M / r + ur^2)
+  #   dρdur = - C1 / (r^2 * ur^2)
+  #   dhdρ = K * Γ * ρ^(Γ - 2)
+  #   du0dur = ur / u0
+  #   return 2 * (u0^2 * h * dhdρ * dρdur + h^2 * u0 * du0dur)
+  # end
 
   # radial domain over which we solve
   rs = collect(LinRange(rmin, rmax, N_rs))
@@ -65,52 +72,275 @@ function solve()
     # - use the previous result for the remaining radi
     # The following converges to the correct solution:
     if r > 2*M
-      ur_guess = sqrt(1 - 2 * M / r) # asymptotic value for ur
+      ur_guess = urc_sign*sqrt(1 - 2 * M / r) # = sqrt(1+gtt), =^= asymptotic value for ur when u0 -> 1
     else
       ur_guess = urs[idx+1] # use previous iterate where asymptotic value is not defined
     end
     F_at_r(ur) = F(ur, r)
-    dF_at_r(ur) = dF(ur, r)
+    # dF_at_r(ur) = dF(ur, r)
     ur = try
-      Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton())
+      # Roots.find_zero((F_at_r,dF_at_r), ur_guess, Roots.Newton())
+      Roots.find_zero(F_at_r, ur_guess)
     catch err
       println()
-      error("Could not determine solution for 'ur' at radius r = $r. " *
-            "Remaining error: $err")
+      @info("Could not determine solution for 'ur' at radius r = $r. ")
+      rethrow(err)
     end
     # @show (ur, r)
-    if ur < 0.0
-      @warn "Aborting due to ur becoming negative at r = $r." 
-      @info "$(N_rs - idx) points skipped ..."
-      break
-    end
+    # if ur < 0.0
+    #   @warn "Aborting due to ur becoming negative at r = $r."
+    #   @info "$(N_rs - idx) points skipped ..."
+    #   break
+    # end
     urs[idx] = ur
   end
 
   ρs = @. C1 / (rs^2 * urs)
   ps = @. K * ρs^Γ
   ϵs = @. ps / ((Γ - 1) * ρs)
-  βr = @. 2 * M / (rs + 2 * M)
+  βr = @. 2 * M / (rs + 2 * M) # = β^r
   α = @. sqrt( rs / (rs + 2 * M) )
-
-  gtt_uu = @. - (2 * M + rs) / rs
-  gtr_uu = @. 2 * M / rs
-  grr_uu = @. (rs - 2 * M) / rs
+  γrr   = @. 1 + 2*M/rs
+  γrr_uu = @. 1/γrr
+  gtt_uu = @. -1/α^2
+  gtr_uu = @. βr/α^2
+  grr_uu = @. γrr_uu - βr*βr/α^2
+  gtt   = @. -α^2 + γrr*βr*βr
+  gtr   = @. γrr * βr
+  grr   = @. γrr
   u0s_d = @. sqrt(1 - 2 * M / rs + urs^2)
   u0s_u = @. gtt_uu * u0s_d + gtr_uu / grr_uu * (urs - gtr_uu * u0s_d)
+  # u0s_u_v2 = @. (gtt_uu * u0s_d + gtr_uu*grr*urs)/(1-gtr*gtr_uu)
+  # @show(maximum(abs, u0s_u .- u0s_u_v2))
   vrs = @. urs / (α * u0s_u) + βr / α
+  vrs_d = @. γrr * vrs
+  hs = @. 1 + K * Γ * ρs^(Γ - 1) / (Γ - 1)
+  res_C1 = @. rs^2 * ρs * urs - C1
+  res_C2 = @. rs^2 * ρs * hs * u0s_d * urs - C2
+  res_C3 = @. hs^2 * u0s_d^2 - C3
+  vs = @. sqrt(γrr * vrs^2)
+
+  return (; rs, urs, vrs, ρs, ps, ϵs, res_C1, res_C2, res_C3, vs, vrs_d,
+            rmin, rmax, rc, Vc, ρc, Γ, M, N_rs, K)
+
+end
+
+
+function solve_ref()
+  # reference data from bugner's solver SphericalAccretion.c (extracted from bamps)
+  M = 1.0
+  rmin, rmax = 1.8*M, 200*M
+  rc, ρc, Γ = 200*M, 7e-4, 4/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+function solve_papadopoulus_1998()
+  # arXiv:gr-qc/9803087v1
+  # works with Eddington-Finkelstein coords (same form as Bugner; but Bugner calls them Kerr-Schild?)
+  M = 1.0
+  rmin, rmax = 0.5*M, 50*M
+  rc, ρc, Γ = 400*M, 1e-2, 4/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+function solve_papadopoulus_1999()
+  # arXiv:gr-qc/9902018v1
+  # works with ingoing Eddington-Finkelstein coords, but the equations appear to be the same
+  # compare with the 1998 paper, in these coordinates G(r) = 0
+  M = 1.0
+  rmin, rmax = 1.5*M, 30*M
+  rc, ρc, Γ = 200*M, 1e-2, 5/3
+  solve(; rmin, rmax, rc, ρc, Γ, M)
+end
+
+
+function save(fname, data)
+  # save
+  filename = joinpath(@__DIR__, fname)
+  h5open(filename, "w") do file
+    @info "Saving results to '$filename' ..."
+    g = create_group(file, "bondi_accretion")
+    HDF5.attributes(g)["rc"]   = data.rc
+    HDF5.attributes(g)["ρc"]   = data.ρc
+    HDF5.attributes(g)["Vc"]   = data.Vc
+    HDF5.attributes(g)["Γ"]    = data.Γ
+    HDF5.attributes(g)["M"]    = data.M
+    HDF5.attributes(g)["K"]    = data.K
+    HDF5.attributes(g)["rmin"] = data.rmin
+    HDF5.attributes(g)["rmax"] = data.rmax
+    g["r"]  = data.rs
+    g["ur"] = data.urs
+    g["vr"] = data.vrs
+    g["ρ"]  = data.ρs
+    g["p"]  = data.ps
+    g["ϵ"]  = data.ϵs
+  end
+  return
+end
+
+
+function plot(data)
+
+  @unpack M, ρs, ps, ϵs, vrs, rs, res_C1, res_C2, res_C3, vrs_d = data
 
   fig = Figure()
-  ax_ρ   = Axis(fig[1,1], xlabel="r", ylabel="ρ")
-  ax_p   = Axis(fig[1,2], xlabel="r", ylabel="p")
-  ax_ϵ   = Axis(fig[2,1], xlabel="r", ylabel="ϵ")
-  ax_v_r = Axis(fig[2,2], xlabel="r", ylabel="v_r")
+  ax_ρ      = Axis(fig[1,1], xlabel=L"r", ylabel=L"\rho")
+  ax_p      = Axis(fig[1,2], xlabel=L"r", ylabel=L"p")
+  ax_ϵ      = Axis(fig[2,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_v_r    = Axis(fig[2,2], xlabel=L"r", ylabel=L"v^r")
+  # ax_res_C1 = Axis(fig[3,1], xlabel=L"r", ylabel=L"R_{C_1}")
+  # ax_res_C2 = Axis(fig[3,2], xlabel=L"r", ylabel=L"R_{C_2}")
+  # ax_res_C3 = Axis(fig[4,1], xlabel=L"r", ylabel=L"R_{C_3}")
+  # ax_v_r_d  = Axis(fig[4,2], xlabel=L"r", ylabel=L"v_r")
+
+  rhor = 2*M # radius of horizon
 
   lines!(ax_ρ, rs, ρs)
   lines!(ax_p, rs, ps)
   lines!(ax_ϵ, rs, ϵs)
   lines!(ax_v_r, rs, vrs)
+  # lines!(ax_res_C1, rs, res_C1)
+  # lines!(ax_res_C2, rs, res_C2)
+  # lines!(ax_res_C3, rs, res_C3)
+  # lines!(ax_v_r_d, rs, vrs_d)
+  vlines!(ax_ρ, rhor, linestyle=:dash)
+  vlines!(ax_p, rhor, linestyle=:dash)
+  vlines!(ax_ϵ, rhor, linestyle=:dash)
+  vlines!(ax_v_r, rhor, linestyle=:dash)
+  # vlines!(ax_v_r_d, rhor, linestyle=:dash)
+
+  # xlims!(ax_ρ, (1.8,10.0))
+  # xlims!(ax_p, (1.8,10.0))
+  # xlims!(ax_ϵ, (1.8,10.0))
+  # xlims!(ax_v_r, (1.8,10.0))
 
   # use save("filename.pdf", fig) to save to pdf
   return fig
 end
+
+
+function plot_compare_ref()
+
+  mydata = solve_ref()
+  csv = readdlm(joinpath(@__DIR__, "bondi-acc-ref.txt"))
+  refdata = (; rs=csv[:,1],
+               ρs=csv[:,2],
+               ϵs=csv[:,3],
+               ps=csv[:,4],
+               vrs=csv[:,5],
+               urs=csv[:,6],
+               ut2=csv[:,7])
+
+  fig = Figure()
+  ax_ρ      = Axis(fig[1,1], xlabel=L"r", ylabel=L"\rho")
+  ax_p      = Axis(fig[1,2], xlabel=L"r", ylabel=L"p")
+  ax_ϵ      = Axis(fig[2,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_vr     = Axis(fig[2,2], xlabel=L"r", ylabel=L"v^r")
+  ax_ur     = Axis(fig[3,1], xlabel=L"r", ylabel=L"u^r")
+
+  lines!(ax_ρ,   mydata.rs,  mydata.ρs)
+  lines!(ax_p,   mydata.rs,  mydata.ps)
+  lines!(ax_ϵ,   mydata.rs,  mydata.ϵs)
+  lines!(ax_vr,  mydata.rs,  mydata.vrs)
+  lines!(ax_ur,  mydata.rs,  mydata.urs)
+  lines!(ax_ρ,   refdata.rs, refdata.ρs, linestyle=:dash)
+  lines!(ax_p,   refdata.rs, refdata.ps, linestyle=:dash)
+  lines!(ax_ϵ,   refdata.rs, refdata.ϵs, linestyle=:dash)
+  lines!(ax_vr,  refdata.rs, refdata.vrs, linestyle=:dash)
+  lines!(ax_ur,  refdata.rs, refdata.urs, linestyle=:dash)
+
+  rhor = 2*mydata.M
+  vlines!(ax_ρ, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_p, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_ϵ, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_vr, rhor, linestyle=:dash, color=:grey)
+  vlines!(ax_ur, rhor, linestyle=:dash, color=:grey)
+
+  return fig
+
+end
+
+
+function plot_compare_papadopoulus_1998()
+
+  mydata = solve_papadopoulus_1998()
+  img_p = load("bondi-acc_papadopoulos1998_p.png")
+  img_ρ = load("bondi-acc_papadopoulos1998_rho.png")
+  # v = sqrt(γ_rr) v^r; is called total velocity in the paper (hidden in text)
+  # the label in the plot is just 'Velocity'
+  img_v = load("bondi-acc_papadopoulos1998_v.png")
+
+  xmin, xmax = 0.0, 15.0
+  ymin_ρ, ymax_ρ = 0, 120.0
+  ymin_p, ymax_p = 0, 2.5
+  ymin_vr, ymax_vr = -0.3, -0.2
+
+  fig = Figure()
+  ax_p  = Axis(fig[1,1], xlabel=L"r", ylabel=L"p")
+  ax_ρ  = Axis(fig[2,1], xlabel=L"r", ylabel=L"\rho")
+  ax_vr = Axis(fig[3,1], xlabel=L"r", ylabel=L"v^r")
+
+  image!(ax_ρ, xmin..xmax, ymin_ρ..ymax_ρ, rotr90(img_ρ))
+  image!(ax_p, xmin..xmax, ymin_p..ymax_p, rotr90(img_p))
+  image!(ax_vr, xmin..xmax, ymin_vr..ymax_vr, rotr90(img_v))
+
+  lines!(ax_ρ,  mydata.rs, mydata.ρs, linestyle=:dash, color=:orange)
+  lines!(ax_p,  mydata.rs, mydata.ps, linestyle=:dash, color=:orange)
+  # don't understand why they plot with a minus here
+  # there seem to be ingoing and outgoing solutions (urc_sign parameter),
+  # but the value for v never fits
+  lines!(ax_vr, mydata.rs, -mydata.vs, linestyle=:dash, color=:orange)
+
+  xlims!(ax_ρ,   (xmin,xmax))
+  xlims!(ax_p,   (xmin,xmax))
+  xlims!(ax_vr,  (xmin,xmax))
+  ylims!(ax_ρ,   (ymin_ρ,ymax_ρ))
+  ylims!(ax_p,   (ymin_p,ymax_p))
+  ylims!(ax_vr,  (ymin_vr,ymax_vr))
+
+  return fig
+end
+
+
+function plot_compare_papadopoulus_1999()
+
+  mydata = solve_papadopoulus_1999()
+  img_ϵ = load("bondi-acc_papadopoulos1999_eps.png")
+  img_ρ = load("bondi-acc_papadopoulos1999_rho.png")
+  img_vr = load("bondi-acc_papadopoulos1999_vr.png")
+
+  xmin, xmax = 1.0, 10.0
+  ymin_ρ, ymax_ρ = 0, 10.0
+  ymin_ϵ, ymax_ϵ = 0, 0.25
+  ymin_vr, ymax_vr = -1.0, 0.0
+
+  fig = Figure()
+  ax_ϵ  = Axis(fig[1,1], xlabel=L"r", ylabel=L"\epsilon")
+  ax_ρ  = Axis(fig[2,1], xlabel=L"r", ylabel=L"\rho")
+  ax_vr = Axis(fig[3,1], xlabel=L"r", ylabel=L"u^r")
+
+  image!(ax_ρ, xmin..xmax, ymin_ρ..ymax_ρ, rotr90(img_ρ))
+  image!(ax_ϵ, xmin..xmax, ymin_ϵ..ymax_ϵ, rotr90(img_ϵ))
+  image!(ax_vr, xmin..xmax, ymin_vr..ymax_vr, rotr90(img_vr))
+
+  lines!(ax_ρ,  mydata.rs, mydata.ρs, linestyle=:dash, color=:orange)
+  lines!(ax_ϵ,  mydata.rs, mydata.ϵs, linestyle=:dash, color=:orange)
+  # don't understand why they plot with a minus here
+  # there seem to be ingoing and outgoing solutions (urc_sign parameter),
+  # but the value for ur never fits
+  lines!(ax_vr, mydata.rs, -mydata.urs, linestyle=:dash, color=:orange)
+  display(extrema(mydata.urs))
+
+  xlims!(ax_ρ,   (xmin,xmax))
+  xlims!(ax_ϵ,   (xmin,xmax))
+  xlims!(ax_vr,  (xmin,xmax))
+  ylims!(ax_ρ,   (ymin_ρ,ymax_ρ))
+  ylims!(ax_ϵ,   (ymin_ϵ,ymax_ϵ))
+  ylims!(ax_vr,  (ymin_vr,ymax_vr))
+
+  return fig
+end
+
+
+
+end # module

initialdata/TOV.jl

@@ -2,58 +2,65 @@ export TOV
 module TOV
 
 
-using dg1d
-using LaTeXStrings
-using Plots
+import CairoMakie
+import CairoMakie: @L_str
+using HDF5
+using OrdinaryDiffEq
 using Printf
 
 
 # array indices of state vector
-const idx_p, idx_m, idx_ϕ, idx_R, idx_r = 1, 2, 3, 4, 5
+const idx_p, idx_m, idx_ϕ, idx_λ = 1, 2, 3, 4
 
 
-function TOV_rhs!(du, u, r, eos::Polytrope)
+function TOV_rhs!(du, u, p, r)
+
+  # parameters for polytropic equation of state
+  (; K, Γ) = p
 
   if r == 0.0
     du .= 0.0
   else
     p, m = u[idx_p], u[idx_m]
 
-    ϵ = try
-      rho = eos(Density, Pressure, p)
-      eps = eos(InternalEnergy, Density, rho)
-      eos(TotalEnergy, Density, rho, InternalEnergy, eps)
+    e = try
+      ρ = (p/K)^(1/Γ)
+      ϵ = K / (Γ-1) * ρ^(Γ-1)
+      ρ*(1+ϵ)
     catch err
       if err isa DomainError
         0.0
       else
-        rethrow(e)
+        rethrow(err)
       end
     end
 
     # dp/dr
-    du[idx_p] = - (ϵ + p) * (m + 4.0 * π * r^3 * p) / (r * (r - 2.0 * m))
+    du[idx_p] = - (e + p) * (m + 4.0 * π * r^3 * p) / (r * (r - 2.0 * m))
     # dm/dr
-    du[idx_m] = 4.0 * π * r^2 * ϵ
+    du[idx_m] = 4.0 * π * r^2 * e
     # dϕ/dr
-    du[idx_ϕ] = - du[idx_p] / (ϵ + p)
-    # dR/dr
-    du[idx_R] = (1.0 - sqrt(1.0 - 2.0 * m / r)) / sqrt(r^2 - 2.0 * m * r)
+    du[idx_ϕ] = - du[idx_p] / (e + p)
+    # regular term of dR/dr ODE
+    du[idx_λ] = (1.0 - sqrt(1.0 - 2.0 * m / r)) / sqrt(r^2 - 2.0 * m * r)
   end
 
   return
 end
 
 
-function solve()
+"""
+Compute TOV solution for a polytropic equation of state.
+  - Γ: polytropic exponent
+  - K: polytropic constant
+  - ρc: central density
+  - rmax: [optional] maximum integration radius, increase when integration does not proceed till surface
+"""
+function solve(; Γ, K, ρc, rmax=30.0)
 
-  # equation of state
-  Γ = 2
-  K = 100
-  # ρc = 8.00 * 1e-3 # central rest mass density
-  ρc = 1.28 * 1e-3 # central rest mass density
-  eos = Polytrope(K = K, Gamma = Γ)
-  pc = eos(Pressure, Density, ρc, InternalEnergy, nothing)
+  # equation of state parameters
+  prms = (; K, Γ)
+  pc = K*ρc^Γ
 
   # initial values
   # Note: We put ϕ(0) = 0 as initial condition for the potential in the interior
@@ -70,155 +77,112 @@ function solve()
   u0[idx_p] = pc
   u0[idx_m] = 0.0
   u0[idx_ϕ] = 0.0
-  u0[idx_R] = 0.0
-
+  u0[idx_λ] = 0.0
 
   dr = 1e-2 # inital step width for r
-  r_range = [0.0, 30.0] # max integration range
-  ODE!(du, u, r) = TOV_rhs!(du, u, r, eos)
-
-
-  # post hook setup to 
-  #   - store solution
-  #   - terminate ODE solver if pressure drops below p_threshold
-  # solution = [u0..., r_range[1]]
-  storage_increment = 1000
-  n_variables = 5
-  solution = zeros(Float64, storage_increment, n_variables)
-  solution[1,:] .= [ u0..., r_range[1] ]
-  counter = 1
-  p_threshold = 1e-14
-  function post_hook(u, r)
-    counter += 1
-    if counter >= size(solution)[1]
-      zs = zeros(Float64, storage_increment, n_variables)
-      solution = vcat(solution, zs)
-    end
-    solution[counter, :] .= [ u..., r ]
-    p = u[idx_p]
-    return (p > p_threshold) # continue as long as this is true
+  r_range = [0.0, rmax] # max integration range
+  prob = ODEProblem(TOV_rhs!,u0,r_range,prms)
+
+  sol = OrdinaryDiffEq.solve(prob, Tsit5(), dt=dr, adaptive=false)
+
+  nsteps = length(sol)
+  p, m, ϕ, λ, r = Float64[], Float64[], Float64[], Float64[], Float64[]
+  for (i,(ui,ri)) in enumerate(zip(sol.u,sol.t))
+    ui[idx_p] < 0 && break
+    push!(p, ui[idx_p])
+    push!(m, ui[idx_m])
+    push!(ϕ, ui[idx_ϕ])
+    push!(λ, ui[idx_λ])
+    push!(r, ri)
   end
 
-  # solve ODE
-  pblm = dg1d.ODEProblem(:rk4)
-  dg1d.evolve_adaptive_ERK(pblm, ODE!, u0, dr, r_range[2], post_hook=post_hook)
+  ρ = @. (p/K)^(1/Γ)
+  ϵ = @. K / (Γ-1) * ρ^(Γ-1)
 
-  # chop off trailing zeros
-  solution = solution[1:counter,:]
+  # attach Schwarzschild solution exterior to star's surface
+  M      = m[end]
+  r_surf = r[end]
+  λ_surf = λ[end]
+  # potential difference (wrt Schwarzschild solution) at star's surface
+  Δϕ_surf = ϕ[end] - (0.5 * log(1.0 - 2.0 * M / r_surf))
+  # adjust interior potential to align with exterior at star's surface
+  ϕ .-= Δϕ_surf
 
-  M      = solution[end, idx_m]
-  r_surf = solution[end, idx_r]
-  R_surf = solution[end, idx_R]
-  @info "M = $M"
-  @info "r = $r_surf"
-  @info "R = $R_surf"
+  R = @. 0.5*r/r_surf * (sqrt(r_surf^2 - 2*M*r_surf) + r_surf - M) * exp(λ - λ_surf)
+  R_surf = R[end]
+  @show M, r_surf, R_surf
+
+  return (; K, Γ, ρc, M, dr, p, m, ϕ, R, ρ, ϵ, r)
+end
 
-  # potential at star's surface
-  Δϕ_surf = solution[end,idx_ϕ] - (0.5 * log(1.0 - 2.0 * M / r_surf))
-  # adjust interior potential to align with exterior at star's surface
-  solution[:,idx_ϕ] .-= Δϕ_surf
 
-  # fix integration constant between isotropic and schwarzschild coordinate
-  Rs = R_surf
-  rs = r_surf
-  r = solution[:,idx_r]
-  R = solution[:,idx_R]
-  @. solution[:,idx_R] = 0.5 * r / rs * (sqrt(rs^2 - 2 * M * rs) + rs - M) * exp(R - Rs)
+function solve_bugner_stable()
+  solve(Γ=2.0, K=100.0, ρc=1.28e-3)
+end
+function solve_bugner_unstable()
+  solve(Γ=2.0, K=100.0, ρc=8e-3)
+end
 
+
+function save(fname, data)
+  @assert endswith(fname, ".h5")
   # save
-  output_dir = dg1d.get_output_dir()
-  filename = joinpath(output_dir, "TOV.h5")
+  filename = joinpath(@__DIR__, fname)
   h5open(filename, "w") do file
     @info "Saving results to '$filename' ..."
     g = create_group(file, "TOV")
-    attributes(g)["Γ"]  = Γ
-    attributes(g)["K"]  = K
-    attributes(g)["ρc"] = ρc
-    attributes(g)["M"]  = M
-    g["ϕ"] = solution[:,idx_ϕ]
-    g["p"] = solution[:,idx_p]
-    g["m"] = solution[:,idx_m]
-    g["R"] = solution[:,idx_R]
-    g["r"] = solution[:,idx_r]
+    attributes(g)["Γ"]    = data.Γ
+    attributes(g)["K"]    = data.K
+    attributes(g)["ρc"]   = data.ρc
+    attributes(g)["M"]    = data.M
+    attributes(g)["dr"]   = data.dr
+    g["ϕ"] = data.ϕ
+    g["p"] = data.p
+    g["m"] = data.m
+    g["R"] = data.R
+    g["r"] = data.r
+    g["ρ"] = data.ρ
+    g["ϵ"] = data.ϵ
   end
-
   return
 end
 
 
+function plot(data)
 
-function read_data(filename)
-  return h5open(filename, "r") do file
-    g = file["TOV"]
-    Γ = read_attribute(g, "Γ")
-    K = read_attribute(g, "K")
-    ρc = read_attribute(g, "ρc")
-    M = read_attribute(g, "M")
-    ϕ = read(g, "ϕ")
-    p = read(g, "p")
-    m = read(g, "m")
-    R = read(g, "R")
-    r = read(g, "r")
-    (; Γ, K, ρc, M, ϕ, p, m, R, r)
-  end
-end
+  (; ρc, dr, M, r, R, ρ, ϵ, m, p, ϕ) = data
+
+  r_surf = r[end]
+  R_surf = R[end]
+  rmax = r_surf*3
+
+  ϕint = ϕ
+  Λint = @. -0.5 * log(1.0-2.0*m/r)
+
+  rext = collect(range(r[end], rmax, step=dr))
+  ϕext = @. 0.5 * log(1.0-2.0*M/rext)
+  Λext = -ϕext
+
+  fig = CairoMakie.Figure(title="TOV star, $(@sprintf "ρc = %.5f, M = %.5f, R=%.5f" ρc M R_surf)")
+  ax_metric = CairoMakie.Axis(fig[1,1], xlabel=L"r")
+  ax_eos    = CairoMakie.Axis(fig[2,1], xlabel=L"r")
+  CairoMakie.linkxaxes!(ax_metric, ax_eos)
 
+  CairoMakie.lines!(ax_metric, r, @.(exp(2.0*ϕint)), label=L"e^{2\phi_{int}}")
+  CairoMakie.lines!(ax_metric, r, @.(exp(-2.0*Λint)), label=L"e^{2\Lambda_{int}}")
+  CairoMakie.lines!(ax_metric, rext, @.(exp(2.0*ϕext)), label=L"e^{2\phi_{ext}}")
+  CairoMakie.lines!(ax_metric, rext, @.(exp(-2.0*Λext)), label=L"e^{2\Lambda_{ext}}")
+  CairoMakie.vlines!(ax_metric, r_surf, linestyle=:dash, color=:black, label=L"r_{surf}")
 
-function plot_potentials()
-
-  output_dir = dg1d.get_output_dir()
-  filename = joinpath(output_dir, "TOV.h5")
-  Γ, K, ρc, M, ϕ, p, m, R, r = read_data(filename)
-
-  eos = Polytrope(K = K, Gamma = Γ)
-  ρ   = @. eos(Density, InternalEnergy, nothing, Pressure, p)
-  ϵ   = @. eos(InternalEnergy, Density, ρ, Pressure, p)
-
-  r_ext = collect(range(r[end], 15.0, length=100))
-  M = m[end]
-  ϕ_ext = @. 0.5 * log(1.0 - 2.0 * M / r_ext)
-
-  ϕ_int = ϕ
-  Λ_int = @. - 0.5 * log(1.0 - 2.0 * m / r)
-  Λ_ext = - ϕ_ext
-
-  Plots.gr()
-  plt_metric = Plots.plot()
-  # Plots.plot(r, 100.0 * ρ, label="100 * ρ(r)")
-  # Plots.plot!(plt, r, 100.0 * ϵ, label="100 * ϵ(r)")
-  # Plots.plot!(plt, r, 100.0 * p, label="100 * p(r)")
-  # Plots.plot!(plt, r, ϕ_int, label="ϕ_int(r)")
-  # Plots.plot!(plt, r_ext, ϕ_ext, label="ϕ_ext(r)")
-  # Plots.plot!(plt, r, Λ_int, label="Λ_int(r)")
-  # Plots.plot!(plt, r_ext, Λ_ext, label="Λ_ext(r)")
-  # Plots.plot!(plt, r, @.(r/R), label="R")
-  # R_ext = @. 1/2 * (sqrt(r_ext^2 - 2 * M * r_ext) + r_ext - M)
-  # Plots.plot!(plt, r_ext, R_ext, label="R")
-  # Plots.plot!(plt, r, m, label="m")
-
-  #### potentials are continuous and differential at surface
-  Plots.plot!(plt_metric, r, @.( exp(2.0*ϕ_int) ), label=L"e^{2 \Phi_{int}}")
-  Plots.plot!(plt_metric, r, @.( exp(-2.0*Λ_int) ), label=L"e^{2 \Lambda_{int}}")
-  Plots.plot!(plt_metric, r_ext, @.( exp(2.0*ϕ_ext) ), label=L"e^{2 \Phi_{ext}}")
-  Plots.plot!(plt_metric, r_ext, @.( exp(-2.0*Λ_ext) ), label=L"e^{2 \Lambda_{ext}}")
-  rs = r[end]
-  Rs = R[end]
-  Plots.plot!(plt_metric, [rs, rs], [0,1], line=(:dash, :black), label=L"r=r_S")
-
-  plt_eos = Plots.plot()
-  Plots.plot!(plt_eos, r, p, label=L"p")
-  Plots.plot!(plt_eos, r, ρ, label=L"\rho")
-  # Plots.plot!(plt_eos, r, ϵ, label="ϵ")
-
-  plt = Plots.plot(plt_metric, plt_eos, layout=(2,1),
-                   xlabel="r",
-                   plot_title="TOV star, $(@sprintf "ρc = %.5f, M = %.5f, R = %.5f" ρc M Rs)",
-                   legendfontsize = 20, guidefontsize = 20,
-                   tickfontsize = 20, titlefontsize = 20,
-                   size=(1200,1200), xlims=(0.0, 15.0))
-  Plots.gui(plt)
+  CairoMakie.lines!(ax_eos, r, p.*100, label=L"100 p")
+  CairoMakie.lines!(ax_eos, r, ρ.*10, label=L"10 \rho")
+  CairoMakie.lines!(ax_eos, r, ϵ, label=L"\epsilon")
+  CairoMakie.vlines!(ax_eos, r_surf, linestyle=:dash, color=:black, label=L"r_{surf}")
 
+  CairoMakie.axislegend(ax_metric, orientation=:horizontal, position=:rb)
+  CairoMakie.axislegend(ax_eos, orientation=:horizontal, position=:rb)
 
+  return fig
 end