From 25ea5e796954471e91c277bff6e75013b7a00c33 Mon Sep 17 00:00:00 2001
From: Nathan Musoke <nathan.musoke@unh.edu>
Date: Wed, 21 Jul 2021 17:01:51 -0300
Subject: [PATCH] feat: phase difference options

Allow user to choose the maximum phase difference and density threshold.
---
 Project.toml      |  2 +-
 src/UltraDark.jl  | 11 ++++-------
 src/config.jl     | 18 +++++++++++++++++-
 src/phase_diff.jl | 25 ++++++++++++++++++++-----
 4 files changed, 42 insertions(+), 14 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3baeaa1..399965e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "UltraDark"
 uuid = "1c8d022d-dfc0-4b41-80ab-3fc7e88cdfea"
 authors = ["Nathan Musoke <n.musoke@auckland.ac.nz>"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
diff --git a/src/UltraDark.jl b/src/UltraDark.jl
index 1fd6d16..8dc2f28 100644
--- a/src/UltraDark.jl
+++ b/src/UltraDark.jl
@@ -12,8 +12,6 @@ export Grids, PencilGrids
 export Config, SimulationConfig, constant_scale_factor
 export OutputConfig
 
-const PHASE_GRAD_LIMIT = π / 4
-
 include("grids.jl")
 include("pencil_grids.jl")
 include("phase_diff.jl")
@@ -151,7 +149,7 @@ function evolve_to!(t_start, t_end, grids, output_config, config::Config.Simulat
 
     while (t < t_end) && ~(t ≈ t_end)
 
-        if max_normed_phase_grad(grids) > PHASE_GRAD_LIMIT
+        if ~good_phase_diff(grids, config)
             output_grids(grids, output_config, -1)
             throw("Phase gradient is too large to continue")
         end
@@ -184,9 +182,8 @@ function simulate(grids, options::Config.SimulationConfig, output_config::Output
 
     # Output initial conditions
     output_grids(grids, output_config, 1)
-
-    if max_normed_phase_grad(grids) > PHASE_GRAD_LIMIT
-        throw("Phase gradient is too large to start")
+    if ~good_phase_diff(grids, options)
+        throw("Initial phase gradient is already in untrusted regime")
     end
 
     for (index, t_end) in enumerate(output_config.output_times[2:end])
@@ -202,7 +199,7 @@ function simulate(grids, options::Config.SimulationConfig, output_config::Output
         output_grids(grids, output_config, index + 1)
     end
 
-    if max_normed_phase_grad(grids) > PHASE_GRAD_LIMIT
+    if ~good_phase_diff(grids, options)
         throw("Phase gradient is too large to end")
     end
 
diff --git a/src/config.jl b/src/config.jl
index 2416878..295fd5d 100644
--- a/src/config.jl
+++ b/src/config.jl
@@ -6,6 +6,8 @@ struct SimulationConfig
     time_step_update_period::Int64
     "function defining a(t)"
     a
+    phase_grad_limit::Float64
+    density_threshold::Float64
 end
 
 function constant_scale_factor(t)
@@ -13,7 +15,21 @@ function constant_scale_factor(t)
 end
 
 function SimulationConfig(time_step_update_period)
-    SimulationConfig(time_step_update_period, constant_scale_factor)
+    SimulationConfig(time_step_update_period, constant_scale_factor, π/4, 1e-6)
+end
+
+function SimulationConfig(time_step_update_period, a)
+    SimulationConfig(time_step_update_period, a, π/4, 1e-6)
+end
+
+function SimulationConfig(
+                          ;
+                          time_step_update_period=10,
+                          a=constant_scale_factor,
+                          phase_grad_limit=π/4,
+                          density_threshold=1e-6
+                         )
+    SimulationConfig(time_step_update_period, constant_scale_factor, phase_grad_limit, density_threshold)
 end
 
 end # module
diff --git a/src/phase_diff.jl b/src/phase_diff.jl
index 3881955..21f9685 100644
--- a/src/phase_diff.jl
+++ b/src/phase_diff.jl
@@ -40,9 +40,8 @@ Compute maximum phase gradient of a grid
 Normalised to ignore large gradients in regions with low density.  These tend
 to be anomalous.
 """
-function max_normed_phase_grad(grids)
-    DENSITY_THRESHOLD = 1e-6
-    tmp = similar(grids.ψx, Float64)
+function max_normed_phase_diff(psi, rho, density_threshold)
+    tmp = similar(psi, Float64)
     n = ndims(tmp)
 
     max_grads = zeros(n)
@@ -52,11 +51,27 @@ function max_normed_phase_grad(grids)
 
     for i in 1:n
         dir = circshift(shift, i)
-        tmp .= abs.(phase_diff(grids.ψx, dir))
-        tmp[grids.ρx / maximum(grids.ρx) .< DENSITY_THRESHOLD] .= NaN
+        tmp .= abs.(phase_diff(psi, dir))
+        tmp[rho / maximum(rho) .< density_threshold] .= NaN
         max_grads[i] = maximum(tmp)
     end
 
     maximum(max_grads)
 
 end
+
+"""
+    bad_phase_diff(grids, config)
+
+Check if the phase of the ψ field is in a trustable regime.
+
+Returns
+"""
+function good_phase_diff(grids, config)
+    density_threshold = config.density_threshold
+    if max_normed_phase_diff(grids.ψx, grids.ρx, density_threshold) > config.phase_grad_limit
+        false
+    else
+        true
+    end
+end