Phase tracing options (#360)

* Added RK45 options to tracePhase

* Options just in freeEnergy

* Tweak

* Added phaseTracerFirstStep option to config

* Fixed typo. Small revert.

* Typo

* Removed unused import

* Added phaseTracerFirstStep to config reader

Models/InertDoubletModel/inertDoubletModelConfig.ini

@@ -80,6 +80,9 @@ tmax = 1.2
 # Desired accuracy of the phase tracer and the resulting FreeEnergy interpolation.
 phaseTracerTol = 1e-8
 
+# First step size in units of the maximum step size. Use None for default algorithm.
+phaseTracerFirstStep = None
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility

Models/ManySinglets/manySingletsConfig.ini

@@ -80,6 +80,9 @@ tmax = 1.2
 # Desired accuracy of the phase tracer and the resulting FreeEnergy interpolation.
 phaseTracerTol = 1e-6
 
+# First step size in units of the maximum step size. Use None for default algorithm.
+phaseTracerFirstStep = None
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility

Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini

@@ -80,6 +80,9 @@ tmax = 1.2
 # Desired accuracy of the phase tracer and the resulting FreeEnergy interpolation.
 phaseTracerTol = 1e-6
 
+# First step size in units of the maximum step size. Use None for default algorithm.
+phaseTracerFirstStep = None
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility

Models/StandardModel/standardModelConfig.ini

@@ -80,6 +80,9 @@ tmax = 1.2
 # Desired accuracy of the phase tracer and the resulting FreeEnergy interpolation.
 phaseTracerTol = 1e-6
 
+# First step size in units of the maximum step size. Use None for default algorithm.
+phaseTracerFirstStep = None
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility

src/WallGo/config.py

@@ -112,6 +112,13 @@ class ConfigThermodynamics:
     Desired accuracy of the phase tracer and the resulting FreeEnergy interpolation.
     """
 
+    phaseTracerFirstStep: float | None = None
+    r"""
+    Starting step for phaseTrace. If a float, this gives the starting step
+    size in units of the maximum step size :py:data:`dT`. If :py:data:`None` then
+    uses the initial step size algorithm of :py:mod:`scipy.integrate.solve_ivp`.
+    """
+
 @dataclass
 class ConfigBoltzmannSolver:
     """ Holds the config of the BoltzmannSolver class. """
@@ -262,7 +269,19 @@ def loadConfigFromFile(self, filePath: str) -> None:
             if 'phaseTracerTol' in keys:
                 self.configThermodynamics.phaseTracerTol = parser.getfloat(
                     "Thermodynamics",
-                    "phaseTracerTol")
+                    "phaseTracerTol"
+                )
+            if 'phaseTracerFirstStep' in keys:
+                # either float or None
+                try:
+                    self.configThermodynamics.phaseTracerFirstStep = parser.getfloat(
+                        "Thermodynamics", "phaseTracerFirstStep"
+                    )
+                except ValueError as valErr:
+                    if str(valErr) == "could not convert string to float: 'None'":
+                        self.configThermodynamics.phaseTracerFirstStep = None
+                    else:
+                        raise
 
         # Read the BoltzmannSolver configs
         if 'BoltzmannSolver' in parser.sections():

src/WallGo/freeEnergy.py

@@ -2,7 +2,6 @@
 Class that does phase tracing, computes the effective potential in the minimum and
 interpolate it.
 """
-
 from dataclasses import dataclass
 import logging
 import numpy as np
@@ -246,6 +245,7 @@ def tracePhase(
         rTol: float = 1e-6,
         spinodal: bool = True,  # Stop tracing if a mass squared turns negative
         paranoid: bool = True,  # Re-solve minimum after every step
+        phaseTracerFirstStep: float | None = None,  # Starting step
     ) -> None:
         r"""Traces minimum of potential
 
@@ -266,12 +266,13 @@ def tracePhase(
         dT : float
             Maximal temperature step size used by the phase tracer.
         rTol : float, optional
-            Relative tolerance of the phase tracing. The default is 1e-6.
+            Relative tolerance of the phase tracing. The default is :py:const:`1e-6`.
         spinodal : bool, optional
             If True, stop tracing if a mass squared turns negative. The default is True.
         paranoid : bool, optional
             If True, re-solve minimum after every step. The default is True.
-
+        phaseTracerFirstStep : float or None, optional
+            If a float, this gives the starting step size in units of the maximum step size :py:data:`dT`. If :py:data:`None` then uses the initial step size algorithm of :py:mod:`scipy.integrate.solve_ivp`. Default is :py:data:`None`
         """
         # make sure the initial conditions are extra accurate
         extraTol = 0.01 * rTol
@@ -325,6 +326,14 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float:
         TMin = max(self.minPossibleTemperature[0], TMin)
         TMax = min(self.maxPossibleTemperature[0], TMax)
 
+        # kwargs for scipy.integrate.solve_ivp
+        scipyKwargs = {
+            "rtol": rTol,
+            "atol": tolAbsolute,
+            "max_step": dT,
+            "first_step": phaseTracerFirstStep,
+        }
+
         # iterating over up and down integration directions
         endpoints = [TMax, TMin]
         for direction in [0, 1]:
@@ -334,10 +343,7 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float:
                 T0,
                 phase0,
                 TEnd,
-                rtol=rTol,
-                atol=tolAbsolute,
-                max_step=dT,
-                first_step=dT,
+                **scipyKwargs,
             )
             while ode.status == "running":
                 try:

src/WallGo/manager.py

@@ -364,8 +364,20 @@ def initTemperatureRange(self) -> None:
         fHighT = self.thermodynamics.freeEnergyHigh
         fLowT = self.thermodynamics.freeEnergyLow
 
-        fHighT.tracePhase(TMinHighT, TMaxHighT, dT, phaseTracerTol)
-        fLowT.tracePhase(TMinLowT, TMaxLowT, dT, phaseTracerTol)
+        fHighT.tracePhase(
+            TMinHighT,
+            TMaxHighT,
+            dT,
+            rTol=phaseTracerTol,
+            phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep,
+        )
+        fLowT.tracePhase(
+            TMinLowT,
+            TMaxLowT,
+            dT,
+            rTol=phaseTracerTol,
+            phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep,
+        )
 
     def setPathToCollisionData(self, directoryPath: pathlib.Path) -> None:
         """