"""
Defines the WallGoManager class which initializes the different object needed for the
wall velocity calculation.
"""
from dataclasses import dataclass
import pathlib
import logging
import numpy as np
# WallGo imports
import WallGo
from .boltzmann import BoltzmannSolver, ETruncationOption
from .containers import PhaseInfo
from .equationOfMotion import EOM
from .exceptions import WallGoError, WallGoPhaseValidationError
from .genericModel import GenericModel
from .grid3Scales import Grid3Scales
from .hydrodynamics import Hydrodynamics
from .hydrodynamicsTemplateModel import HydrodynamicsTemplateModel
from .thermodynamics import Thermodynamics
from .results import WallGoResults
from .config import Config
[docs]
@dataclass
class WallSolverSettings:
"""
Settings for the WallSolver.
"""
bIncludeOffEquilibrium: bool = True
"""If False, will ignore all out-of-equilibrium effects (no Boltzmann solving).
"""
meanFreePathScale: float = 50.0
"""Estimate of the mean free path of the plasma in :math:`1/T_n`. This will be used
to set the tail lengths in the Grid object. Default is 100.
"""
wallThicknessGuess: float = 5.0
r"""
Initial guess of the wall thickness that will be used to solve the EOM, in units
:math:`1/T_n`. Default is 5.
"""
[docs]
@dataclass
class WallSolver:
"""
Data class containing classes and settings for the wall velocity computation.
"""
eom: EOM
"""EOM object"""
grid: Grid3Scales
"""Grid3Scales object used to describe the mapping between the physical
coordinates to the compact ones."""
boltzmannSolver: BoltzmannSolver
""" BoltzmannSolver object used to solve the Boltzmann equation."""
initialWallThickness: float
"""Initial wall thickness used by the solver. Should be expressed in physical
units (the units used in :py:class:`WallGo.EffectivePotential`)."""
[docs]
class WallGoManager:
"""Manages WallGo program flow
The WallGoManager is a 'control' class which collects together and manages
all the various parts of the WallGo Python package for the computation of
the bubble wall velocity.
"""
[docs]
def __init__(self) -> None:
""""""
# Initialise the configs with the default values
self.config = Config()
# Set the default verbosity level to logging.INFO
self.setVerbosity(logging.INFO)
# default to working directory
self.collisionDirectory = pathlib.Path.cwd()
# These we currently have to keep cached, otherwise we can't construct
# a sensible WallSolver:
## TODO init these to None or have other easy way of checking if they
## have been properly initialized
self.model: GenericModel
self.hydrodynamics: Hydrodynamics
self.phasesAtTn: PhaseInfo
self.thermodynamics: Thermodynamics
[docs]
def getMomentumGridSize(self) -> int:
"""
Returns the momentum grid size.
"""
return self.config.configGrid.momentumGridSize
[docs]
def setVerbosity(self, verbosityLevel: int) -> None:
"""
Set the verbosity level.
Parameters
----------
verbosityLevel : int
Verbosity level. Follows the standard convention of the logging module where
:py:const:`DEBUG=10`, :py:const:`INFO=20`, :py:const:`WARNING=30` and :py:const:`ERROR=40`. In WallGo, most of the
information is shown at the :py:const:`INFO` level. At the :py:const:`DEBUG` level, more
information about the calculation of the pressure at each iteration is
shown.
"""
logging.basicConfig(format="%(message)s", level=verbosityLevel, force=True)
[docs]
def setupThermodynamicsHydrodynamics(
self,
phaseInfo: WallGo.PhaseInfo,
veffDerivativeScales: WallGo.VeffDerivativeSettings,
freeEnergyArraysHighT: WallGo.FreeEnergyArrays = None,
freeEnergyArraysLowT: WallGo.FreeEnergyArrays = None,
) -> None:
r"""Must run before :py:meth:`solveWall()` and companions.
Initialization of internal objects related to equilibrium thermodynamics and
hydrodynamics. Specifically, we verify that the input PhaseInfo is valid
(distinct phases can be found in the effective potential),
estimate the relevant temperature range for wall solving and create efficient
approximations for phase free energies over this range using interpolation.
Finally, it initializes :py:class:`Hydrodynamics` and confirms that it can find a
reasonable value for the Jouguet velocity (the transition from hybdrid to
detonation solutions).
You are required to run this function whenever details of your
physics model change to keep the manager's internal state up to date.
Parameters
----------
phaseInfo : PhaseInfo
WallGo object containing the approximate positions of the minima and
nucleation temperature.
veffDerativeScales : VeffDerivativeSettings
WallGo dataclass containing the typical temprature and field scale over
which the potential varies.
Returns
-------
"""
assert (
phaseInfo.phaseLocation1.numFields() == self.model.fieldCount
and phaseInfo.phaseLocation2.numFields() == self.model.fieldCount
), "Invalid PhaseInfo input, field counts don't match those defined in model"
self.model.getEffectivePotential().configureDerivatives(veffDerivativeScales)
# Checks that phase input makes sense with the user-specified Veff
self.validatePhaseInput(phaseInfo)
self.initTemperatureRange(
freeEnergyArraysHighT=freeEnergyArraysHighT,
freeEnergyArraysLowT=freeEnergyArraysLowT,
)
## Should we write these to a result struct?
logging.info("Temperature ranges:")
logging.info(
"High-T phase: TMin = "
f"{self.thermodynamics.freeEnergyHigh.minPossibleTemperature[0]}, "
f"TMax = {self.thermodynamics.freeEnergyHigh.maxPossibleTemperature[0]}"
)
logging.info(
"Low-T phase: TMin = "
f"{self.thermodynamics.freeEnergyLow.minPossibleTemperature[0]}, "
f"TMax = {self.thermodynamics.freeEnergyLow.maxPossibleTemperature[0]}"
)
self.thermodynamics.setExtrapolate()
self._initHydrodynamics(self.thermodynamics)
if (
not np.isfinite(self.hydrodynamics.vJ)
or self.hydrodynamics.vJ > 1
or self.hydrodynamics.vJ < 0
):
raise WallGoError(
"Failed to solve Jouguet velocity at input temperature!",
data={
"vJ": self.hydrodynamics.vJ,
"temperature": phaseInfo.temperature,
},
)
logging.info(f"Jouguet: {self.hydrodynamics.vJ}")
# TODO return some results struct
[docs]
def isModelValid(self) -> bool:
"""True if a valid model is currently registered."""
return self.model is not None
[docs]
def registerModel(self, model: GenericModel) -> None:
"""
Register a physics model with WallGo.
Parameters
----------
model : GenericModel
GenericModel object that describes the model studied.
"""
assert isinstance(model, GenericModel)
assert (
model.fieldCount > 0
), "WallGo model must contain at least one classical field"
self.model = model
[docs]
def initTemperatureRange(
self,
freeEnergyArraysHighT: WallGo.FreeEnergyArrays = None,
freeEnergyArraysLowT: WallGo.FreeEnergyArrays = None,
) -> None:
r"""
Determine the relevant temperature range and trace the phases
over this range. Interpolate the free energy in both phases and
store in internal :py:class:`Thermodynamics` object.
Parameters
----------
freeEnergyArraysHighT : WallGo.FreeEnergyArrays, optional
If provided, use these arrays to initialize the high-T free energy object.
If None, the phase will be traced.
freeEnergyArraysLowT : WallGo.FreeEnergyArrays, optional
If provided, use these arrays to initialize the low-T free energy object.
If None, the phase will be traced.
"""
assert self.phasesAtTn is not None
assert self.isModelValid()
assert (
self.model.getEffectivePotential().areDerivativesConfigured()
), "Must have called effectivePotential.configureDerivatives()"
Tn = self.phasesAtTn.temperature
self.thermodynamics = Thermodynamics(
self.model.getEffectivePotential(),
Tn,
self.phasesAtTn.phaseLocation2,
self.phasesAtTn.phaseLocation1,
)
# Let's turn these off so that things are more transparent
self.thermodynamics.freeEnergyHigh.disableAdaptiveInterpolation()
self.thermodynamics.freeEnergyLow.disableAdaptiveInterpolation()
try:
# Use the template model to find an estimate of the minimum and maximum
# required temperature. We do not solve hydrodynamics inside the bubble, so
# we are only interested in T- (the temperature right at the wall).
hydrodynamicsTemplate = HydrodynamicsTemplateModel(self.thermodynamics)
logging.info(
f"vwLTE in the template model: {hydrodynamicsTemplate.findvwLTE()}"
)
except WallGoError as error:
# Throw new error with more info
raise WallGoPhaseValidationError(
error.message, self.phasesAtTn, error.data
) from error
# Raise an error if this is an inverse PT (if epsilon is negative)
if hydrodynamicsTemplate.epsilon < 0:
raise WallGoError(
f"WallGo requires epsilon={hydrodynamicsTemplate.epsilon} to be "
"positive."
)
phaseTracerTol = self.config.configThermodynamics.phaseTracerTol
# Estimate of the dT needed to reach the desired tolerance considering
# the error of a cubic spline scales like dT**4.
dT = (
self.model.getEffectivePotential().derivativeSettings.temperatureVariationScale
* phaseTracerTol**0.25
)
# Construct high and low temperature free energy objects
fHighT = self.thermodynamics.freeEnergyHigh
fLowT = self.thermodynamics.freeEnergyLow
# Try to construct interpolations if arrays are given
loadedHigh = False
loadedLow = False
if freeEnergyArraysHighT is not None:
# If the user provided free energy arrays, use them to initialize the
# free energy objects.
try:
fHighT.constructInterpolationFromArray(
freeEnergyArraysHighT,
dT,
)
loadedHigh = True
logging.info("Using user-provided high-T free energy arrays.")
except (ValueError, WallGoError) as e:
raise WallGoError(
f"Failed to load high-T free energy arrays: \n {e}"
) from e
if freeEnergyArraysLowT is not None:
# If the user provided free energy arrays, use them to initialize the
# free energy objects.
try:
fLowT.constructInterpolationFromArray(
freeEnergyArraysLowT,
dT,
)
loadedLow = True
logging.info("Using user-provided low-T free energy arrays.")
except (ValueError, WallGoError) as e:
raise WallGoError(
f"Failed to load high-T free energy arrays: \n {e}"
) from e
# If the user did not provide free energy arrays, we trace the phases
if loadedHigh and loadedLow:
return
# Maximum values for T+ and T- are reached at the Jouguet velocity
_, _, THighTMaxTemplate, TLowTMaxTemplate = hydrodynamicsTemplate.findMatching(
0.99 * hydrodynamicsTemplate.vJ
)
# Minimum value for T- is reached at small wall velocity. The minimum value
# for T+ is the nucleation temperature.
_, _, TLowTMinTemplate, _ = hydrodynamicsTemplate.findMatching(1e-3)
if THighTMaxTemplate is None:
THighTMaxTemplate = self.config.configHydrodynamics.tmax * Tn
if TLowTMaxTemplate is None:
TLowTMaxTemplate = self.config.configHydrodynamics.tmax * Tn
if TLowTMinTemplate is None:
TLowTMinTemplate = self.config.configHydrodynamics.tmin * Tn
phaseTracerTol = self.config.configThermodynamics.phaseTracerTol
interpolationDegree = self.config.configThermodynamics.interpolationDegree
# Estimate of the dT needed to reach the desired tolerance considering
# the error of a cubic spline scales like dT**4.
dT = (
self.model.getEffectivePotential().derivativeSettings.temperatureVariationScale
* phaseTracerTol**0.25
)
"""Since the template model is an approximation of the full model,
and since the temperature profile in the wall could be non-monotonous,
we should not take exactly the TMin and TMax from the template model.
We use a configuration parameter to determine the TMin and TMax that we
use in the phase tracing.
"""
TMinHighT = Tn * self.config.configThermodynamics.tmin
TMaxHighT = THighTMaxTemplate * self.config.configThermodynamics.tmax
TMinLowT = TLowTMinTemplate * self.config.configThermodynamics.tmin
TMaxLowT = TLowTMaxTemplate * self.config.configThermodynamics.tmax
# Only trace if the corresponding file wasn't loaded
if not loadedHigh:
fHighT.tracePhase(
TMinHighT,
TMaxHighT,
dT,
rTol=phaseTracerTol,
phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep,
)
if not loadedLow:
fLowT.tracePhase(
TMinLowT,
TMaxLowT,
dT,
rTol=phaseTracerTol,
phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep,
)
[docs]
def setPathToCollisionData(self, directoryPath: pathlib.Path) -> None:
"""
Specify path to collision files for use with the Boltzmann solver.
This does not necessarily load the files immediately.
Args:
directoryPath (pathlib.Path): Directory containing the .hdf5 collision data.
Returns:
None
"""
# TODO validate? or not if we want to allow the data to be generated on the fly
# should at least validate that it is not an existing file
self.collisionDirectory = directoryPath
[docs]
def getCurrentCollisionDirectory(self) -> pathlib.Path:
"""
Returns the path to the directory with the collision files.
"""
return self.collisionDirectory
[docs]
def wallSpeedLTE(self) -> float:
"""
Solves wall speed in the Local Thermal Equilibrium (LTE) approximation.
Returns
-------
float
Wall velocity in LTE.
"""
return self.hydrodynamics.findvwLTE()
[docs]
def solveWall(
self,
wallSolverSettings: WallSolverSettings,
) -> WallGoResults:
r"""
Solves for the wall velocity
Solves the coupled scalar equation of motion and the Boltzmann equation.
Must be ran after :py:meth:`analyzeHydrodynamics()` because
the solver depends on thermodynamical and hydrodynamical
info stored internally in the :py:class:`WallGoManager`.
Parameters
----------
wallSolverSettings : WallSolverSettings
Configuration settings for the solver.
Returns
-------
WallGoResults
Object containing the wall velocity and EOM solution, as well as different
quantities used to assess the accuracy of the solution.
"""
solver: WallSolver = self.setupWallSolver(wallSolverSettings)
return solver.eom.findWallVelocityDeflagrationHybrid(
solver.initialWallThickness
)
[docs]
def solveWallDetonation(
self,
wallSolverSettings: WallSolverSettings,
onlySmallest: bool = True,
) -> list[WallGoResults]:
"""
Finds all the detonation solutions by computing the pressure on a grid
and interpolating to find the roots.
Parameters
----------
wallSolverSettings : WallSolverSettings
Configuration settings for the solver.
onlySmallest : bool, optional
Whether or not to only look for one solution. If True, the solver will
stop the calculation after finding the first root. If False, it will
continue looking for solutions until it reaches the maximal velocity.
Returns
-------
list[WallGoResults]
List containing the detonation solutions. If no solutions were found,
returns a wall velocity of 0 if the pressure is always positive, or 1 if
it is negative (runaway wall). If it is positive at vmin and negative at
vmax, the outcome is uncertain and would require a time-dependent analysis,
so it returns an empty list.
"""
solver: WallSolver = self.setupWallSolver(wallSolverSettings)
assert solver.initialWallThickness
rtol = self.config.configEOM.errTol
nbrPointsMin = self.config.configEOM.nbrPointsMinDeton
nbrPointsMax = self.config.configEOM.nbrPointsMaxDeton
overshootProb = self.config.configEOM.overshootProbDeton
vmin = max(self.hydrodynamics.vJ + 1e-3, self.hydrodynamics.slowestDeton())
vmax = self.config.configEOM.vwMaxDeton
if vmin >= vmax:
raise WallGoError(
"In WallGoManager.solveWallDetonation(): vmax must be larger than vmin",
{"vmin": vmin, "vmax": vmax},
)
return solver.eom.findWallVelocityDetonation(
vmin,
vmax,
solver.initialWallThickness,
nbrPointsMin,
nbrPointsMax,
overshootProb,
rtol,
onlySmallest,
)
[docs]
def setupWallSolver(self, wallSolverSettings: WallSolverSettings) -> WallSolver:
r"""Helper for constructing an :py:class:`EOM` object whose state is tied to equilibrium
and hydrodynamical information stored in the :py:class:`WallGoManager`.
Specifically, uses results of :py:meth:`setupThermodynamicsHydrodynamics()` to find
optimal grid settings for wall solving, and creates :py:class:`Grid` and :py:class:`BoltzmannSolver`
objects to use from within the :py:class:`EOM`. Be aware that the created :py:class:`EOM` object can
change state if its creator :py:class:`WallGoManager` instance is modified
(e.g. if :py:meth:`setupThermodynamicsHydrodynamics()` is called)!
"""
assert (
self.phasesAtTn.temperature is not None
and self.isModelValid()
and self.hydrodynamics is not None
), "Run WallGoManager.setupThermodynamicsHydrodynamics() before wall solving"
Tnucl: float = self.phasesAtTn.temperature
gridMomentumFalloffScale = Tnucl
wallThickness = wallSolverSettings.wallThicknessGuess
meanFreePathScale = wallSolverSettings.meanFreePathScale
grid: Grid3Scales = self.buildGrid(
wallThickness,
meanFreePathScale,
gridMomentumFalloffScale,
)
# Factor that multiplies the collision term in the Boltzmann equation.
collisionMultiplier = self.config.configBoltzmannSolver.collisionMultiplier
truncationOption = ETruncationOption[
self.config.configBoltzmannSolver.truncationOption
]
# Hardcode basis types here: Cardinal for z, Chebyshev for pz, pp
boltzmannSolver = BoltzmannSolver(
grid,
basisM="Cardinal",
basisN="Chebyshev",
collisionMultiplier=collisionMultiplier,
truncationOption=truncationOption,
)
boltzmannSolver.updateParticleList(self.model.outOfEquilibriumParticles)
bShouldLoadCollisions = wallSolverSettings.bIncludeOffEquilibrium
if bShouldLoadCollisions:
# This throws if collision load fails, let caller handle the exception.
# TODO may be cleaner to handle it here and return an invalid solver
boltzmannSolver.loadCollisions(self.collisionDirectory)
eom: EOM = self.buildEOM(grid, boltzmannSolver, meanFreePathScale)
eom.includeOffEq = wallSolverSettings.bIncludeOffEquilibrium
return WallSolver(eom, grid, boltzmannSolver, wallThickness / Tnucl)
def _initHydrodynamics(self, thermodynamics: Thermodynamics) -> None:
r"""
Initialize the :py:class:`Hydrodynamics` object.
Parameters
----------
thermodynamics : Thermodynamics
Thermodynamics object.
"""
tmax = self.config.configHydrodynamics.tmax
tmin = self.config.configHydrodynamics.tmin
rtol = self.config.configHydrodynamics.relativeTol
atol = self.config.configHydrodynamics.absoluteTol
self.hydrodynamics = Hydrodynamics(thermodynamics, tmax, tmin, rtol, atol)
[docs]
def buildGrid(
self,
wallThicknessIni: float,
meanFreePathScale: float,
initialMomentumFalloffScale: float,
) -> Grid3Scales:
r"""
Initialize a :py:class:`Grid3Scales` object
Parameters
----------
wallThicknessIni : float
Initial guess of the wall thickness that will be used to solve the EOM.
Should be expressed in units of :math:`1/T_n`.
meanFreePathScale : float
Estimate of the mean free path of the plasma. This will be used to set the
tail lengths in the Grid object. Should be expressed in units of :math:`1/T_n`.
initialMomentumFalloffScale : float
TODO documentation. Should be close to temperature at the wall
"""
gridN = self.config.configGrid.momentumGridSize
gridM = self.config.configGrid.spatialGridSize
ratioPointsWall = self.config.configGrid.ratioPointsWall
smoothing = self.config.configGrid.smoothing
Tnucl = self.phasesAtTn.temperature
# We divide by Tnucl to get it in physical units of length
tailLength = (
max(
meanFreePathScale,
0.5 * wallThicknessIni * (1.0 + 3.0 * smoothing) / ratioPointsWall,
)
/ Tnucl
)
if gridN % 2 == 0:
raise ValueError(
"You have chosen an even number N of momentum-grid points. "
"WallGo only works with odd N, please change it to an odd number."
)
return Grid3Scales(
gridM,
gridN,
tailLength,
tailLength,
wallThicknessIni / Tnucl,
initialMomentumFalloffScale,
ratioPointsWall,
smoothing,
)
[docs]
def buildEOM(
self,
grid: Grid3Scales,
boltzmannSolver: BoltzmannSolver,
meanFreePathScale: float,
) -> EOM:
r"""
Constructs an :py:class:`EOM` object using internal state from the :py:class:`WallGoManager`,
and the input :py:class:`Grid` and :py:class:`BoltzmannSolver`.
Parameters
----------
grid : Grid3Scales
Grid3Scales object used to describe the mapping between the physical
coordinates to the compact ones.
boltzmannSolver : BoltzmannSolver
BoltzmannSolver object used to solve the Boltzmann equation.
meanFreePathScale : float
Estimate of the mean free path of the plasma. This will be used to set the
tail lengths in the Grid object. Should be expressed in :math:`1/T_n`.
Returns
-------
EOM
EOM object.
"""
numberOfFields = self.model.fieldCount
errTol = self.config.configEOM.errTol
maxIterations = self.config.configEOM.maxIterations
pressRelErrTol = self.config.configEOM.pressRelErrTol
conserveEnergy = self.config.configEOM.conserveEnergyMomentum
wallThicknessBounds = self.config.configEOM.wallThicknessBounds
wallOffsetBounds = self.config.configEOM.wallOffsetBounds
Tnucl = self.phasesAtTn.temperature
return EOM(
boltzmannSolver,
self.thermodynamics,
self.hydrodynamics,
grid,
numberOfFields,
meanFreePathScale / Tnucl, # We divide by Tnucl to get physical units
wallThicknessBounds,
wallOffsetBounds,
includeOffEq=True,
forceEnergyConservation=conserveEnergy,
forceImproveConvergence=False,
errTol=errTol,
maxIterations=maxIterations,
pressRelErrTol=pressRelErrTol,
)