added r(rstar, M)

Manifest.toml

@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.3"
 manifest_format = "2.0"
-project_hash = "b10dce6173fef906c5a8327d99ba626a58dbf940"
+project_hash = "b123ba64d6d16309c98d31f425d679e6c86a913c"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
@@ -83,7 +83,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.1+0"
 
 [[deps.CpuId]]
 deps = ["Markdown"]
@@ -183,6 +183,11 @@ git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
 version = "1.5.0"
 
+[[deps.LambertW]]
+git-tree-sha1 = "c5ffc834de5d61d00d2b0e18c96267cffc21f648"
+uuid = "984bce1d-4616-540c-a9ee-88d1112d94c9"
+version = "0.4.6"
+
 [[deps.LayoutPointers]]
 deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"]
 git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277"
@@ -331,7 +336,7 @@ weakdeps = ["Adapt"]
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]

Project.toml

@@ -7,6 +7,7 @@ version = "1.0.0-DEV"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridFunctions = "28d702cb-9102-4da8-a82a-5023545b1c0c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 

src/InitialData.jl

@@ -1,6 +1,15 @@
 module InitialData
 
 using ForwardDiff
+using LambertW
+
+function rstar(r::T, M::T) where {T}
+    return r + 2M * log(abs(r / (2M) - one(T)))
+end
+
+function r(rstar::T, M::T) where {T}
+    return 2(M + M * lambertw(exp(-one(T) + rstar / (2M))))
+end
 
 @inline function Gaussian1D(t::Real, r::Real, σ::Real, r0::Real)
     return exp(-(r - r0 + t)^2 / σ^2)
@@ -15,16 +24,15 @@ end
 end
 
 
-
 #creating general functions
 @inline function Lamb(l::T, r::T, M::T) where {T}
     return (l - one(T)) * (l + 2one(T)) + 6M / r
 end
 @inline function f(r::T, M::T) where {T}
     return one(T) - 2M / r
 end
-@inline function mul(l::Real)
-    return (l - 1) * (l - 2)
+@inline function mul(l::T) where {T}
+    return (l - one(T)) * (l - 2one(T))
 end
 
 #defining the function
@@ -121,6 +129,7 @@ end
 end
 
 # Potential for the function
+# THIS IS EVALUATED at r and not rstar
 @inline function Vpot(l::Real, r::Real, M::Real)
     return f(r, M) / (Lamb(l, r, M)^2) * (mul(l)^2 * ((mul(l) + 2) / r^2 + 6 * M / r^3) + 36 * M^2 / r^4 * (mul(l) + 2 * M / r))
 end

src/ODE.jl

@@ -77,9 +77,8 @@ using ..InitialData
                            $boundary_type")
                 end
             else
-                ## you need to provide the transformation from rstar to r
-                ## to get the correct potential
-                dtΠ[i] = drΨ + InitialData.Vpot(l, rstar[i], M) * Φ[i]
+                ri = InitialData.r(rstar[i], M) # transform rstar to r for potential
+                dtΠ[i] = drΨ + InitialData.Vpot(l, ri, M) * Φ[i]
                 dtΨ[i] = drΠ
             end
         end

src/Run.jl

@@ -7,24 +7,16 @@ using GridFunctions
 
 export zerilli
 
-
-function rstar(r, M)
-    return r + 2M * log(abs(r - 2M))
-end
-
-# you need to provide the transformations
-function r(rstar, M) end
-
-function zerilli(ncells::Integer, tf::Real, cfl::Real;
-    boundary_type::Symbol=:radiative, folder="", save_every=1)
+function zerilli(domain::Vector, ncells::Integer,
+                 tf::Real, cfl::Real, M::Real=1.0, l::Real=0;
+                 boundary_type::Symbol=:radiative,
+                 folder="", save_every=1)
 
     @assert boundary_type === :radiative || boundary_type === :reflective
-    M = 1.0
-    horizon = 2M
-    grid = UniformGrid([horizon, 20M + horizon], ncells)
 
-    N = ncells + 1
+    grid = UniformGrid(domain, ncells)
 
+    N = ncells + 1
     nt = ceil(Int64, tf / (cfl * spacing(grid)))
     t = UniformGrid([0.0, tf], nt)
 
@@ -35,10 +27,9 @@ function zerilli(ncells::Integer, tf::Real, cfl::Real;
     Π = InitialData.dtGaussian1D.(0, coords(grid), σ, r0)
     Ψ = InitialData.drGaussian1D.(0, coords(grid), σ, r0)
 
-    l = 0.0
     params = (h=spacing(grid), N=N, bc=boundary_type, t=t, ti=coords(t),
-        dt=spacing(t), save_every=save_every, M=M, grid=grid,
-        folder=folder, cfl=cfl, rcoord=coords(grid), l=l)
+              dt=spacing(t), save_every=save_every, M=M, grid=grid,
+              folder=folder, cfl=cfl, rcoord=coords(grid), l=l)
     statevector = hcat(Φ, Π, Ψ)
     Integrator.solve(ODE.rhs!, statevector, params)
     return nothing