diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 759d45bc..4c0653ca 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -2,6 +2,20 @@ Changelog ================== +1.1.0 (2025-08-20) +================== + +* Phase tracing updates to improve algorithm stability. +* User may now pass the free energy of a phase as arrays. +* Spectral truncation options added to avoid aliasing, with an automatic method used by default. +* Error estimate for the wall speed updated to make it more theoretically sound. +* Linearization criterion updated so that `linearizationCriterion2` estimates the second order correction. +* Error estimate due to the Tanh ansatz now computed and added to `WallGoResults`. +* Accuracy of energy-momentum conservation computed, and printed for logging level `DEBUG`. +* New tests added for the Standard Model with light Higgs. +* `PTTools `_ example file added. + + 1.0.0 (2024-11-07) ================== diff --git a/Models/InertDoubletModel/inertDoubletModel.py b/Models/InertDoubletModel/inertDoubletModel.py index 1e0c1995..65a963fd 100644 --- a/Models/InertDoubletModel/inertDoubletModel.py +++ b/Models/InertDoubletModel/inertDoubletModel.py @@ -378,7 +378,7 @@ def jCW( # Note that we are taking the absolute value of the mass in the log here, # instead of using EImaginaryOption = ABS_ARGUMENT, because we do not - # want the absolute value in the product of massSq ans rgScale + # want the absolute value in the product of massSq and rgScale return degreesOfFreedom*np.array( massSq * massSq * (np.log(np.abs(massSq / rgScale**2)) - c) + 2 * massSq * rgScale**2 diff --git a/Models/InertDoubletModel/inertDoubletModelConfig.ini b/Models/InertDoubletModel/inertDoubletModelConfig.ini index 0f1bfa96..f517ef40 100644 --- a/Models/InertDoubletModel/inertDoubletModelConfig.ini +++ b/Models/InertDoubletModel/inertDoubletModelConfig.ini @@ -83,6 +83,9 @@ phaseTracerTol = 1e-8 # First step size in units of the maximum step size. Use None for default algorithm. phaseTracerFirstStep = None +# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives. +interpolationDegree = 1 + [BoltzmannSolver] # Factor multiplying the collision term in the Boltzmann equation. # Can be used for testing or for studying the solution's sensibility @@ -96,3 +99,6 @@ basisM = Cardinal # The momentum polynomial basis type, either Cardinal or Chebyshev. basisN = Chebyshev + +# Whether or not to truncate the spectral expansion +truncationOption = AUTO \ No newline at end of file diff --git a/Models/ManySinglets/manySingletsConfig.ini b/Models/ManySinglets/manySingletsConfig.ini index 16065aa6..4a87096d 100644 --- a/Models/ManySinglets/manySingletsConfig.ini +++ b/Models/ManySinglets/manySingletsConfig.ini @@ -83,6 +83,9 @@ phaseTracerTol = 1e-6 # First step size in units of the maximum step size. Use None for default algorithm. phaseTracerFirstStep = None +# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives. +interpolationDegree = 1 + [BoltzmannSolver] # Factor multiplying the collision term in the Boltzmann equation. # Can be used for testing or for studying the solution's sensibility @@ -96,3 +99,6 @@ basisM = Cardinal # The momentum polynomial basis type, either Cardinal or Chebyshev. basisN = Chebyshev + +# Whether or not to truncate the spectral expansion +truncationOption = AUTO \ No newline at end of file diff --git a/Models/SingletStandardModel_Z2/exampleCollisionDefs.py b/Models/SingletStandardModel_Z2/exampleCollisionDefs.py index 18837342..fc2df125 100644 --- a/Models/SingletStandardModel_Z2/exampleCollisionDefs.py +++ b/Models/SingletStandardModel_Z2/exampleCollisionDefs.py @@ -1,3 +1,10 @@ +""" +This Python script, exampleCollisionDefs.py, +uses the implementation of the minimal Standard Model +extension in singletStandardModelZ2.py and adds example +definitions for the Boltzmann collision integrals for +WallGoCollision. +""" import WallGoCollision import pathlib import sys diff --git a/Models/SingletStandardModel_Z2/exampleOutputThermodynamics.py b/Models/SingletStandardModel_Z2/exampleOutputThermodynamics.py new file mode 100644 index 00000000..6ccf9bb6 --- /dev/null +++ b/Models/SingletStandardModel_Z2/exampleOutputThermodynamics.py @@ -0,0 +1,129 @@ +""" +This Python script, exampleOutputThermodynamics.py, +uses the implementation of the minimal Standard Model +extension in singletStandardModelZ2.py and gives methods +for saving the thermodynamics of the model for later use +in e.g. PTTools. +""" + +import sys +import pathlib +import numpy as np +import h5py + +# WallGo imports +import WallGo # Whole package, in particular we get WallGo._initializeInternal() + +# Add the SingletStandardModelZ2 folder to the path to import SingletSMZ2 +sys.path.append(str(pathlib.Path(__file__).resolve().parent)) +from singletStandardModelZ2 import SingletSMZ2 + + +def main() -> None: + + manager = WallGo.WallGoManager() + + # Model definition is done in the SingletSMZ2 class + model = SingletSMZ2() + manager.registerModel(model) + + inputParameters = { + "RGScale": 125.0, + "v0": 246.0, + "MW": 80.379, + "MZ": 91.1876, + "Mt": 173.0, + "g3": 1.2279920495357861, + "mh1": 125.0, + "mh2": 120.0, + "lHS": 0.9, + "lSS": 1.0, + } + + model.updateModel(inputParameters) + + # Creates a Thermodynamics object, containing all thermodynamic functions + # by tracing the high- and low-temperature phases. + # For temperatures outside of the range of existence of the phases, the + # thermodynamic functions are extrapolated using the template model. + manager.setupThermodynamicsHydrodynamics( + WallGo.PhaseInfo( + temperature=100.0, # nucleation temperature + phaseLocation1=WallGo.Fields([0.0, 200.0]), + phaseLocation2=WallGo.Fields([246.0, 0.0]), + ), + WallGo.VeffDerivativeSettings( + temperatureVariationScale=10.0, fieldValueVariationScale=[10.0, 10.0] + ), + ) + + Tcrit = manager.thermodynamics.findCriticalTemperature(dT = 0.1) + + # Generate tables for the thermodynamics functions of both phases. + # Note that if the minimum or maximum temperature chosen is outside of + # the range of existence of the phases, the extrapolation is used. + + # Temperature range and step for high-temperature phase + highTPhaseRange = (80.0, 140.0, 0.1) + # Temperature range and step for low-temperature phase + lowTPhaseRange = (80.0, 110.0, 0.1) + + # Create temperature arrays + temp_high = np.arange(highTPhaseRange[0], highTPhaseRange[1], highTPhaseRange[2]) + temp_low = np.arange(lowTPhaseRange[0], lowTPhaseRange[1], lowTPhaseRange[2]) + + # Evaluate thermodynamic functions for high-temperature phase + p_high = np.array([manager.thermodynamics.pHighT(T) for T in temp_high]) + e_high = np.array([manager.thermodynamics.eHighT(T) for T in temp_high]) + csq_high = np.array([manager.thermodynamics.csqHighT(T) for T in temp_high]) + + # Evaluate thermodynamic functions for low-temperature phase + p_low = np.array([manager.thermodynamics.pLowT(T) for T in temp_low]) + e_low = np.array([manager.thermodynamics.eLowT(T) for T in temp_low]) + csq_low = np.array([manager.thermodynamics.csqLowT(T) for T in temp_low]) + + # Get temperature limits; outside of these limits, extrapolation has been used + max_temp_high = manager.thermodynamics.freeEnergyHigh.maxPossibleTemperature[0] + min_temp_high = manager.thermodynamics.freeEnergyHigh.minPossibleTemperature[0] + max_temp_low = manager.thermodynamics.freeEnergyLow.maxPossibleTemperature[0] + min_temp_low = manager.thermodynamics.freeEnergyLow.minPossibleTemperature[0] + + # Create HDF5 file + filename = pathlib.Path(__file__).resolve().parent/"thermodynamics_data.h5" + + with h5py.File(filename, 'w') as f: + # Global attributes + f.attrs["model_label"] = "singletSMZ2" + f.attrs["critical_temperature"] = Tcrit + f.attrs["nucleation_temperature"] = 100. + + # High-temperature phase group + high_group = f.create_group("high_temperature_phase") + high_group.create_dataset("temperature", data=temp_high) + high_group.create_dataset("pressure", data=p_high) + high_group.create_dataset("energy_density", data=e_high) + high_group.create_dataset("sound_speed_squared", data=csq_high) + high_group.attrs["max_possible_temperature"] = max_temp_high + high_group.attrs["min_possible_temperature"] = min_temp_high + + # Low-temperature phase group + low_group = f.create_group("low_temperature_phase") + low_group.create_dataset("temperature", data=temp_low) + low_group.create_dataset("pressure", data=p_low) + low_group.create_dataset("energy_density", data=e_low) + low_group.create_dataset("sound_speed_squared", data=csq_low) + low_group.attrs["max_possible_temperature"] = max_temp_low + low_group.attrs["min_possible_temperature"] = min_temp_low + + print(f"Thermodynamics data saved to {filename}") + print("Model: singletSMZ2") + print(f"Critical temperature: {Tcrit:.2f} GeV") + print(f"High-T phase: {len(temp_high)} temperature points from {temp_high[0]:.1f} to {temp_high[-1]:.1f} GeV") + print(f"Low-T phase: {len(temp_low)} temperature points from {temp_low[0]:.1f} to {temp_low[-1]:.1f} GeV") + print(f"High-T phase valid range: {min_temp_high:.1f} to {max_temp_high:.1f} GeV") + print(f"Low-T phase valid range: {min_temp_low:.1f} to {max_temp_low:.1f} GeV") + + +## Don't run the main function if imported to another file +if __name__ == "__main__": + main() diff --git a/Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini b/Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini index 8a1bd1c4..57393798 100644 --- a/Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini +++ b/Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini @@ -83,6 +83,9 @@ phaseTracerTol = 1e-6 # First step size in units of the maximum step size. Use None for default algorithm. phaseTracerFirstStep = None +# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives. +interpolationDegree = 1 + [BoltzmannSolver] # Factor multiplying the collision term in the Boltzmann equation. # Can be used for testing or for studying the solution's sensibility @@ -96,3 +99,6 @@ basisM = Cardinal # The momentum polynomial basis type, either Cardinal or Chebyshev. basisN = Chebyshev + +# Whether or not to truncate the spectral expansion +truncationOption = AUTO diff --git a/Models/StandardModel/standardModel.py b/Models/StandardModel/standardModel.py index f9afe3c3..293561a8 100644 --- a/Models/StandardModel/standardModel.py +++ b/Models/StandardModel/standardModel.py @@ -25,7 +25,6 @@ A Study of the electroweak phase transition dynamics, Phys.Rev.D 52 (1995) 7182-7204 doi:10.1103/PhysRevD.52.7182 """ - import sys import pathlib import numpy as np @@ -299,9 +298,9 @@ def evaluate( # pylint: disable=R0914 lambdaT = self.modelParameters["lambda"] - 3 / ( 16 * np.pi * np.pi * self.modelParameters["v0"] ** 4 ) * ( - 2 * mW**4 * np.log(mW**2 / (ab * T**2) + 1e-100) - + mZ**4 * np.log(mZ**2 / (ab * T**2) + 1e-100) - - 4 * mt**4 * np.log(mt**2 / (af * T**2) + 1e-100) + 2 * mW**4 * np.log(mW**2 / (ab * T**2) ) + + mZ**4 * np.log(mZ**2 / (ab * T**2) ) + - 4 * mt**4 * np.log(mt**2 / (af * T**2) ) ) cT: float | np.ndarray = self.modelParameters["C0"] + 1 / ( @@ -317,7 +316,7 @@ def evaluate( # pylint: disable=R0914 potentialT: float | np.ndarray = ( self.modelParameters["D"] * (T**2 - self.modelParameters["T0sq"]) * v**2 - - cT * T**2 * pow(v, 2) * np.log(np.abs(v / T)) + - cT * T**2 * pow(v, 2) * np.log(np.abs(v / T) + 1e-100) # Avoid log(0) - eT * T * pow(v, 3) + lambdaT / 4 * pow(v, 4) ) @@ -529,7 +528,7 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]: of the Higgs mass. """ valuesMH = [0.0, 34.0, 50.0, 70.0, 81.0] - valuesTn = [57.192, 70.579, 83.426, 102.344, 113.575] + valuesTn = [57.1958, 70.5793, 83.4251, 102.344, 113.575] output: list[ExampleInputPoint] = [] @@ -555,8 +554,8 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]: WallGo.WallSolverSettings( # we actually do both cases in the common example bIncludeOffEquilibrium=True, - meanFreePathScale=100.0, # In units of 1/Tnucl - wallThicknessGuess=20.0, # In units of 1/Tnucl + meanFreePathScale=50.0, # In units of 1/Tnucl + wallThicknessGuess=5.0, # In units of 1/Tnucl ), ) ) diff --git a/Models/StandardModel/standardModelConfig.ini b/Models/StandardModel/standardModelConfig.ini index 8a1bd1c4..c129f3b6 100644 --- a/Models/StandardModel/standardModelConfig.ini +++ b/Models/StandardModel/standardModelConfig.ini @@ -16,7 +16,7 @@ momentumGridSize = 11 ratioPointsWall = 0.5 # Smoothing factor of the mapping function (the larger the smoother) -smoothing = 0.1 +smoothing = 0.3 [EquationOfMotion] # The absolute error tolerance for the wall-velocity result @@ -83,6 +83,9 @@ phaseTracerTol = 1e-6 # First step size in units of the maximum step size. Use None for default algorithm. phaseTracerFirstStep = None +# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives. +interpolationDegree = 1 + [BoltzmannSolver] # Factor multiplying the collision term in the Boltzmann equation. # Can be used for testing or for studying the solution's sensibility @@ -96,3 +99,6 @@ basisM = Cardinal # The momentum polynomial basis type, either Cardinal or Chebyshev. basisN = Chebyshev + +# Whether or not to truncate the spectral expansion +truncationOption = AUTO diff --git a/Models/Yukawa/yukawaCollisionForPairs.py b/Models/Yukawa/yukawaCollisionForPairs.py new file mode 100644 index 00000000..afa1417e --- /dev/null +++ b/Models/Yukawa/yukawaCollisionForPairs.py @@ -0,0 +1,164 @@ +""" +Demonstration of the computeIntegralsForPair methods. + +This script shows how to compute collision integrals for specific particle pairs +rather than the full set of out-of-equilibrium partice pairs. +""" + +import sys +import pathlib +import os + +# WallGo imports +import WallGo +import WallGoCollision + +# Add the Yukawa folder to the path to import YukawaModel +sys.path.append(str(pathlib.Path(__file__).resolve().parent)) +from yukawa import YukawaModel + + +def initCollisionModel(wallGoModel: "YukawaModel") -> "WallGoCollision.PhysicsModel": + """Initialize the Collision model and set the seed.""" + + import WallGoCollision # pylint: disable = C0415 + + # Collision integrations utilize Monte Carlo methods, so RNG is involved. + # We can set the global seed for collision integrals as follows. + # This is optional; by default the seed is 0. + WallGoCollision.setSeed(0) + + collisionModelDefinition = ( + WallGo.collisionHelpers.generateCollisionModelDefinition(wallGoModel) + ) + + # Add in-equilibrium particles that appear in collision processes + # The out-of-equilibrium particles are taken from the definition in the model file + phiParticle = WallGoCollision.ParticleDescription() + phiParticle.name = "phi" + phiParticle.index = 0 + phiParticle.bInEquilibrium = True + phiParticle.bUltrarelativistic = True + phiParticle.type = WallGoCollision.EParticleType.eBoson + # mass-sq function not required or used for UR particles, + # and it cannot be field-dependent for collisions. + # Backup of what the vacuum mass was intended to be: + + parameters = WallGoCollision.ModelParameters() + + parameters.add("y", wallGoModel.modelParameters["y"]) + parameters.add("gamma", wallGoModel.modelParameters["gamma"]) + parameters.add("lam", wallGoModel.modelParameters["lam"]) + parameters.add("v", 0.0) + + # fermion asymptotic thermal mass^2 (twice the static thermal mass) + # in units of T + parameters.add( + "mf2", 1 / 8 * wallGoModel.modelParameters["y"] ** 2 + ) + # scalar thermal mass^2 in units of T + parameters.add( + "ms2", + +wallGoModel.modelParameters["lam"] / 24.0 + + wallGoModel.modelParameters["y"] ** 2.0 / 6.0, + ) + + collisionModelDefinition.defineParticleSpecies(phiParticle) + collisionModelDefinition.defineParameters(parameters) + + collisionModel = WallGoCollision.PhysicsModel(collisionModelDefinition) + + return collisionModel + +def demonstrateComputeIntegralsForPair(collision_model: "WallGoCollision.PhysicsModel", save_dir: str | None = None) -> None: + """Demonstrate computing collision integrals for specific particle pairs using computeIntegralsForPair().""" + + gridSize = 3 # Small grid for faster demonstration + + # Set default save directory if not provided + if save_dir is None: + save_dir = f"CollisionOutput_N{gridSize}_Pairs" + + # Create output directory + os.makedirs(save_dir, exist_ok=True) + + # Load matrix elements FIRST (before creating collision tensor) + print("Loading matrix elements...") + matrix_element_file = pathlib.Path(__file__).resolve().parent / "MatrixElements" / "matrixElements.yukawa.json" + bShouldPrintMatrixElements = True + if not collision_model.loadMatrixElements(str(matrix_element_file), bShouldPrintMatrixElements): + print(f"FATAL: Failed to load matrix elements from {matrix_element_file}") + return {} + + # Create collision tensor AFTER loading matrix elements + print("Creating collision tensor...") + collisionTensor = collision_model.createCollisionTensor(gridSize) + + # Configure integration options + options = WallGoCollision.IntegrationOptions() + options.maxIntegrationMomentum = 10 # Reduced for speed + options.absoluteErrorGoal = 1e-6 + options.relativeErrorGoal = 1e-1 + options.maxTries = 20 + options.calls = 10000 # Reduced for speed + options.bIncludeStatisticalErrors = True + + collisionTensor.setIntegrationOptions(options) + + # Set verbosity for demonstration + verbosity = WallGoCollision.CollisionTensorVerbosity() + verbosity.bPrintElapsedTime = True + verbosity.progressReportPercentage = 0.5 + verbosity.bPrintEveryElement = False + collisionTensor.setIntegrationVerbosity(verbosity) + + # Demonstrate computing collision integrals for multiple particle pairs + print("\n5. Computing collision integrals for particle pairs") + + # Define the particle pairs to compute + particle_pairs = [("psiL", "psiL"), ("psiR", "psiR")] + + results = {} + + for particle1, particle2 in particle_pairs: + print(f"\nComputing collision integrals for pair: {particle1} - {particle2}") + + # Use computeIntegralsForPair to compute the collision integrals for this specific pair + pair_result = collisionTensor.computeIntegralsForPair(particle1, particle2) + + if pair_result is not None: + print(f"Successfully computed integrals for {particle1}-{particle2}") + + # Save directly to the main directory + filename = os.path.join(save_dir, f"{particle1}_{particle2}.h5") + print(f"Saving to: {filename}") + + pair_result.writeToHDF5(filename) + print(f"Successfully saved {particle1}-{particle2}") + + results[f"{particle1}_{particle2}"] = pair_result + else: + print(f"Failed to compute integrals for {particle1}-{particle2} (returned None)") + + return + + +if __name__ == "__main__": + + # Initialize the Yukawa model + wallGoModel = YukawaModel() + + # Set example parameters + inputParameters = { + "sigma": 0.0, + "msq": 1.0, + "gamma": -1.2, + "lam": 0.10, + "y": 0.55, + "mf": 0.30, + } + wallGoModel.modelParameters.update(inputParameters) + + collisionModel = initCollisionModel(wallGoModel) + + results = demonstrateComputeIntegralsForPair(collisionModel) diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 2ef751d4..bcc36841 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -10,6 +10,8 @@ Here we collect a number of different example model files. The Yukawa model file examples/wallGoExampleBase examples/standardModel examples/singletScalarExtension + examples/singletScalarExtensionCollisions + examples/singletScalarExtensionThermodynamics examples/manySinglets examples/inertDoubletModel examples/yukawa \ No newline at end of file diff --git a/docs/source/examples/singletScalarExtensionCollisions.rst b/docs/source/examples/singletScalarExtensionCollisions.rst new file mode 100644 index 00000000..05f4f710 --- /dev/null +++ b/docs/source/examples/singletScalarExtensionCollisions.rst @@ -0,0 +1,6 @@ +================================================ +Singlet Scalar Extension: Collisions Definitions +================================================ + +.. literalinclude:: ../../../Models/SingletStandardModel_Z2/exampleCollisionDefs.py + :language: py diff --git a/docs/source/examples/singletScalarExtensionThermodynamics.rst b/docs/source/examples/singletScalarExtensionThermodynamics.rst new file mode 100644 index 00000000..6c763f63 --- /dev/null +++ b/docs/source/examples/singletScalarExtensionThermodynamics.rst @@ -0,0 +1,6 @@ +=============================================== +Singlet Scalar Extension: Thermodynamics Output +=============================================== + +.. literalinclude:: ../../../Models/SingletStandardModel_Z2/exampleOutputThermodynamics.py + :language: py \ No newline at end of file diff --git a/docs/source/faqs.rst b/docs/source/faqs.rst index c2431891..aae71072 100644 --- a/docs/source/faqs.rst +++ b/docs/source/faqs.rst @@ -13,7 +13,7 @@ General - **How should I cite WallGo?** WallGo is free and open source, but if you use WallGo in your work, we ask that you - support us by please citing the WallGo paper, `arXiv:2411.04970 `_. The complete BibTex citation from `Inspire `_ is:: + support us by please citing the WallGo paper, `JHEP 04 (2025) 101 `_. The complete BibTex citation from `Inspire `_ is:: @article{Ekstedt:2024fyq, author = "Ekstedt, Andreas and Gould, Oliver and Hirvonen, Joonas and Laurent, Benoit and Niemi, Lauri and Schicho, Philipp and van de Vis, Jorinde", @@ -22,8 +22,11 @@ General archivePrefix = "arXiv", primaryClass = "hep-ph", reportNumber = "CERN-TH-2024-174, DESY-24-162, HIP-2024-21/TH", - month = "11", - year = "2024" + doi = "10.1007/JHEP04(2025)101", + journal = "JHEP", + volume = "04", + pages = "101", + year = "2025" } @@ -156,6 +159,18 @@ Effective potentials enum :py:class:`WallGo.PotentialTools.EImaginaryOption`. See the docs for more details. +Free energy +----------- +- **I already know the value of the field and the effective potential as a function of temperature, can I provide these to WallGo to circumvent the phase tracing?** + + If the phase tracing does not work properly for your model, or if you want to speed up the + initialization phase, you can provide arrays with the values of the field(s) in the minimum of the + potential and the corresponding effective potential for the appropriate temperature range. + These are passed as a :py:class:`WallGo.FreeEnergyArrays` object, to the function + :py:meth:`WallGo.WallGoManager.setupThermodynamicsHydrodynamics()`. These arrays are optional arguments; + if they are not provided, WallGo will execute its default phase tracing algorithm. + + Settings ======== diff --git a/docs/source/index.rst b/docs/source/index.rst index 196abf8d..1959b7fa 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -2,7 +2,7 @@ WallGo documentation ====================================== -WallGo is an open source code for the computation of the bubble wall velocity and bubble wall width in first-order cosmological phase transitions. If you use WallGo, please cite the WallGo paper `arXiv:2411.04970 `_. :footcite:p:`Ekstedt:2024fyq` +WallGo is an open source code for the computation of the bubble wall velocity and bubble wall width in first-order cosmological phase transitions. If you use WallGo, please cite the WallGo paper `JHEP 04 (2025) 101 `_. :footcite:p:`Ekstedt:2024fyq` As the universe cooled after the Hot Big Bang, it may have gone through any number of cosmological first-order phase transitions. Such transitions proceed via the nucleation and growth of bubbles, as shown in the image below. :footcite:p:`Weir_2016` The collisions of these bubbles may lead to an observable gravitational wave signal today, depending on the speed of the bubble walls as they collide. diff --git a/docs/source/refs.bib b/docs/source/refs.bib index ab915eb4..3e1e1ca6 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -54,7 +54,9 @@ @article{Ekstedt:2024fyq archivePrefix = "arXiv", primaryClass = "hep-ph", reportNumber = "CERN-TH-2024-174, DESY-24-162, HIP-2024-21/TH", - month = "11", - year = "2024", - journal = "arXiv" + doi = "10.1007/JHEP04(2025)101", + journal = "JHEP", + volume = "04", + pages = "101", + year = "2025" } \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index eca889ea..afe3feb4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,6 @@ dependencies = [ "h5py>=3.8.0", "numpy>=1.21.5", "scipy>=1.9.0", - "deprecated>=1.2.14", ] readme = "README.md" diff --git a/src/WallGo/__init__.py b/src/WallGo/__init__.py index dba74acb..3fb522a8 100644 --- a/src/WallGo/__init__.py +++ b/src/WallGo/__init__.py @@ -5,14 +5,14 @@ import importlib # package level modules -from .boltzmann import BoltzmannSolver +from .boltzmann import BoltzmannSolver, ETruncationOption from .config import Config from .collisionArray import CollisionArray -from .containers import PhaseInfo, BoltzmannBackground, BoltzmannDeltas, WallParams +from .containers import PhaseInfo, BoltzmannBackground, BoltzmannDeltas, FreeEnergyArrays, WallParams from .effectivePotential import EffectivePotential, VeffDerivativeSettings from .exceptions import WallGoError, WallGoPhaseValidationError, CollisionLoadError from .fields import Fields -from .freeEnergy import FreeEnergy +from .freeEnergy import FreeEnergy, FreeEnergyValueType from .genericModel import GenericModel from .grid import Grid from .grid3Scales import Grid3Scales @@ -21,10 +21,10 @@ from .interpolatableFunction import InterpolatableFunction, EExtrapolationType from .manager import WallGoManager, WallSolverSettings from .particle import Particle -from .polynomial import Polynomial +from .polynomial import Polynomial, SpectralConvergenceInfo from .thermodynamics import Thermodynamics from .equationOfMotion import EOM -from .results import WallGoResults +from .results import WallGoResults, ESolutionType from .utils import getSafePathToResource # list of submodules for lazy importing diff --git a/src/WallGo/boltzmann.py b/src/WallGo/boltzmann.py index e6668a91..fe1913b7 100644 --- a/src/WallGo/boltzmann.py +++ b/src/WallGo/boltzmann.py @@ -6,11 +6,12 @@ import typing from copy import deepcopy import logging +from enum import Enum, auto import numpy as np import findiff # finite difference methods from .containers import BoltzmannBackground, BoltzmannDeltas from .grid import Grid -from .polynomial import Polynomial +from .polynomial import Polynomial, SpectralConvergenceInfo from .particle import Particle from .collisionArray import CollisionArray from .results import BoltzmannResults @@ -20,6 +21,19 @@ import importlib +class ETruncationOption(Enum): + """Enums for what to do with truncating the spectral expansion.""" + + NONE = auto() + """Do not truncate early, use all coefficients.""" + + AUTO = auto() + """Truncate early if it seems the UV is not converging.""" + + THIRD = auto() + """Drop the last third of the coefficients.""" + + class BoltzmannSolver: """ Class for solving Boltzmann equations for small deviations from equilibrium. @@ -33,6 +47,7 @@ class BoltzmannSolver: offEqParticles: list[Particle] background: BoltzmannBackground collisionArray: CollisionArray + truncationOption: ETruncationOption def __init__( self, @@ -41,6 +56,7 @@ def __init__( basisN: str = "Chebyshev", derivatives: str = "Spectral", collisionMultiplier: float = 1.0, + truncationOption: ETruncationOption = ETruncationOption.AUTO, ): """ Initialisation of BoltzmannSolver @@ -62,6 +78,10 @@ def __init__( collisionMultiplier : float, optional Factor by which the collision term is multiplied. Can be used for testing. Default is 1.0. + truncationOption : ETruncationOption, optional + Option for truncating the spectral expansion. Default is + ETruncationOption.AUTO. Other options + are ETruncationOption.NONE and ETruncationOption.THIRD. Returns ------- @@ -83,8 +103,9 @@ def __init__( self.basisM = basisM # Momentum polynomial type self.basisN = basisN - + self.collisionMultiplier = collisionMultiplier + self.truncationOption = truncationOption # These are set, and can be updated, by our member functions # TODO: are these None types the best way to go? @@ -117,7 +138,10 @@ def updateParticleList(self, offEqParticles: list[Particle]) -> None: self.offEqParticles = offEqParticles - def getDeltas(self, deltaF: typing.Optional[np.ndarray] = None) -> BoltzmannResults: + def getDeltas( + self, + deltaF: typing.Optional[np.ndarray] = None, + ) -> BoltzmannResults: """ Computes Deltas necessary for solving the Higgs equation of motion. @@ -139,11 +163,18 @@ def getDeltas(self, deltaF: typing.Optional[np.ndarray] = None) -> BoltzmannResu if deltaF is None: deltaF = self.solveBoltzmannEquations() - # getting (optimistic) estimate of truncation error - truncationError = self.estimateTruncationError(deltaF) + # checking spectral convergence + deltaF, shapeTruncated, spectralPeaks = self.checkSpectralConvergence(deltaF) - # getting criteria for validity of linearization - criterion1, criterion2 = self.checkLinearization(deltaF) + # getting (optimistic) estimate of truncation error + truncationError = self.estimateTruncationError( + deltaF, shapeTruncated + ) + truncatedTail = ( + shapeTruncated[1] != deltaF.shape[1], + shapeTruncated[2] != deltaF.shape[2], + shapeTruncated[3] != deltaF.shape[3], + ) particles = self.offEqParticles @@ -202,8 +233,8 @@ def getDeltas(self, deltaF: typing.Optional[np.ndarray] = None) -> BoltzmannResu deltaF=deltaF, Deltas=Deltas, truncationError=truncationError, - linearizationCriterion1=criterion1, - linearizationCriterion2=criterion2, + truncatedTail=truncatedTail, + spectralPeaks=spectralPeaks, ) def solveBoltzmannEquations(self) -> np.ndarray: @@ -218,15 +249,15 @@ def solveBoltzmannEquations(self) -> np.ndarray: where :math:`\mathcal{L}` is the Lioville operator, :math:`\mathcal{C}` is the collision operator, and :math:`\mathcal{S}` is the source. - + As regards the indicies, - + - :math:`\alpha, \beta, \gamma` denote points on the coordinate lattice :math:`\{\xi^{(\alpha)},p_{z}^{(\beta)},p_{\Vert}^{(\gamma)}\}`, - :math:`i, j, k` denote elements of the basis of spectral functions :math:`\{\bar{T}_i, \bar{T}_j, \tilde{T}_k\}`, - :math:`a, b` denote particle species. - + For more details see the WallGo paper. Parameters @@ -262,7 +293,7 @@ def solveBoltzmannEquations(self) -> np.ndarray: return deltaF - def estimateTruncationError(self, deltaF: np.ndarray) -> float: + def estimateTruncationError(self, deltaF: np.ndarray, shapeTruncated: tuple[int, ...]) -> float: r""" Quick estimate of the polynomial truncation error using John Boyd's Rule-of-thumb-2: the last coefficient of a Chebyshev @@ -287,38 +318,192 @@ def estimateTruncationError(self, deltaF: np.ndarray) -> float: # sum(|deltaF|) as the norm deltaFPoly.changeBasis(("Array", "Chebyshev", "Chebyshev", "Chebyshev")) - deltaFMeanAbs = np.sum( - np.abs(deltaFPoly.coefficients), + deltaFTuncated = deltaFPoly.coefficients[ + :shapeTruncated[0], + :shapeTruncated[1], + :shapeTruncated[2], + :shapeTruncated[3], + ] + deltaFSumAbs = np.sum( + np.abs(deltaFTuncated), axis=(1, 2, 3), ) # estimating truncation errors in each direction truncationErrorChi = np.sum( - np.abs(deltaFPoly.coefficients[:, -1, :, :]), + np.abs(deltaFTuncated[:, -1, :, :]), axis=(1, 2), - ) + ) / deltaFSumAbs truncationErrorPz = np.sum( - np.abs(deltaFPoly.coefficients[:, :, -1, :]), + np.abs(deltaFTuncated[:, :, -1, :]), axis=(1, 2), - ) + ) / deltaFSumAbs truncationErrorPp = np.sum( - np.abs(deltaFPoly.coefficients[:, :, :, -1]), + np.abs(deltaFTuncated[:, :, :, -1]), axis=(1, 2), - ) + ) / deltaFSumAbs # estimating the total truncation error as the maximum of these three return max( # type: ignore[no-any-return] - np.max(truncationErrorChi / deltaFMeanAbs), - np.max(truncationErrorPz / deltaFMeanAbs), - np.max(truncationErrorPp / deltaFMeanAbs), + np.max(truncationErrorChi), + np.max(truncationErrorPz), + np.max(truncationErrorPp), ) + def checkSpectralConvergence(self, deltaF: np.ndarray) -> tuple[np.ndarray, tuple[int, int, int, int], tuple[int, int, int]]: + """ + Check for spectral convergence. + + Fits to the exponential slope of the last 1/3 of coefficients in the + Chebyshev basis. Also returns the + + Parameters + ---------- + deltaF : array_like + The solution for which to estimate the truncation error, + a rank 3 array, with shape :py:data:`(len(z), len(pz), len(pp))`. + + Returns + ------- + deltaFTruncated : np.ndarray + Potentially truncated version of input :py:data:`deltaF`, padded with zeros if truncated, so same shape as input. + shapeTruncated : tuple[int, int, int, int] + Shape of truncated array. + spectralPeaks : tuple[int, int, int] + Indices of the peaks in the (potentially truncated) spectral expansion. + """ + # constructing Polynomial + basisTypes = ("Array", self.basisM, self.basisN, self.basisN) + basisNames = ("Array", "z", "pz", "pp") + deltaFPoly = Polynomial(deltaF, self.grid, basisTypes, basisNames, False) + truncatedShape = list(deltaF.shape) + + # changing to Chebyshev basis + deltaFPoly.changeBasis(("Array", "Chebyshev", "Chebyshev", "Chebyshev")) + + # looking at convergence of spectral expansion + spectralCoeffsChi = np.sum( + np.abs(deltaFPoly.coefficients), + axis=(0, 2, 3), + ) + spectralCoeffsPz = np.sum( + np.abs(deltaFPoly.coefficients), + axis=(0, 1, 3), + ) + spectralCoeffsPp = np.sum( + np.abs(deltaFPoly.coefficients), + axis=(0, 1, 2), + ) + + # how much to cut, if truncating + cutSpatial = -((self.grid.M - 1) // 3) + cutMomentum = -((self.grid.N - 1) // 3) + + # checking spectral convergence of spatial direction + chiConvergenceTailInfo = SpectralConvergenceInfo( + spectralCoeffsChi[cutSpatial:], + # weightPower=0, + offset=self.grid.M - 1 + cutSpatial, + ) + + # checking spectral convergence of pz direction + pzConvergenceTailInfo = SpectralConvergenceInfo( + spectralCoeffsPz[cutMomentum:], + # weightPower=1, # removed as max(pz) only grows as log(N) + offset=self.grid.N - 1 + cutMomentum, + ) + + # checking spectral convergence of pp direction + ppConvergenceTailInfo = SpectralConvergenceInfo( + spectralCoeffsPp[cutMomentum:], + # weightPower=2, # removed as max(pp) only grows as log(N) + offset=self.grid.N - 1 + cutMomentum, + ) + + allTailsConverging = ( + chiConvergenceTailInfo.apparentConvergence and + pzConvergenceTailInfo.apparentConvergence and + ppConvergenceTailInfo.apparentConvergence + ) + + # Deciding what to do based on truncationOption + if self.truncationOption == ETruncationOption.AUTO: + # if the slope is not definitely negative, we will truncate + if not chiConvergenceTailInfo.apparentConvergence: + deltaFPoly.coefficients[:, cutSpatial:, :, :] = 0 + truncatedShape[1] = deltaF.shape[1] + cutSpatial + if not pzConvergenceTailInfo.apparentConvergence: + deltaFPoly.coefficients[:, :, cutMomentum:, :] = 0 + truncatedShape[2] = deltaF.shape[2] + cutMomentum + if not ppConvergenceTailInfo.apparentConvergence: + deltaFPoly.coefficients[:, :, :, cutMomentum:] = 0 + truncatedShape[3] = deltaF.shape[3] + cutMomentum + elif self.truncationOption == ETruncationOption.THIRD: + # truncating regardless + deltaFPoly.coefficients[:, cutSpatial:, :, :] = 0 + deltaFPoly.coefficients[:, :, cutMomentum:, :] = 0 + deltaFPoly.coefficients[:, :, :, cutMomentum:] = 0 + truncatedShape[1:] = [ + deltaF.shape[1] + cutSpatial, + deltaF.shape[2] + cutMomentum, + deltaF.shape[3] + cutMomentum, + ] + if allTailsConverging: + logging.info( + "Tails of spectral expansions converging but truncated, consider changing truncation option." + ) + else: + # not truncating regardless + if not allTailsConverging: + logging.info( + "Tails of spectral expansions not converging, consider changing truncation option, or changing grid parameters." + ) + + # checking spectral convergence of z direction + chiConvergenceInfo = SpectralConvergenceInfo( + spectralCoeffsChi[:truncatedShape[1]], weightPower=0 + ) + + # checking spectral convergence of pz direction + pzConvergenceInfo = SpectralConvergenceInfo( + spectralCoeffsPz[:truncatedShape[2]], weightPower=1 + ) + + # checking spectral convergence of pp direction + ppConvergenceInfo = SpectralConvergenceInfo( + spectralCoeffsPp[:truncatedShape[3]], weightPower=2 + ) + + # putting together the spectral peaks + spectralPeaks = ( + chiConvergenceInfo.spectralPeak, + pzConvergenceInfo.spectralPeak, + ppConvergenceInfo.spectralPeak, + ) + + if self.truncationOption == ETruncationOption.NONE: + return deltaF, tuple(truncatedShape), spectralPeaks + + # changing back to original basis + deltaFPoly.changeBasis(basisTypes) + + return deltaFPoly.coefficients, tuple(truncatedShape), spectralPeaks + + @staticmethod + def _smoothTruncation(length: int, cut: int, sharp: float = 3) -> np.ndarray: + """ + Internal function to smooth the truncation of the spectral expansion. """ + x = np.arange(length) + return 1 / (1 + np.exp(sharp * (x - cut))) + def checkLinearization( self, deltaF: typing.Optional[np.ndarray] = None - ) -> tuple[np.ndarray, np.ndarray]: + ) -> tuple[float, float]: r""" - Compute two criteria to verify the validity of the linearization of the - Boltzmann equation: :math:`\delta f/f_{eq}` and :math:`C[\delta f]/L[\delta f]`. + Compute two criteria to verify the validity of the linearisation of the + Boltzmann equation: :math:`\delta f/f_{eq}` and + :math:`\delta f_2/(f_{eq}+\delta f)`, with :math:`\delta f_2` the first-order + correction due to nonlinearities. To be valid, at least one of the two criteria must be small for each particle. Parameters @@ -348,20 +533,51 @@ def checkLinearization( ) deltaFPoly.changeBasis(("Array", "Cardinal", "Cardinal", "Cardinal")) + # Computing \delta f^2 + deltaFSqPoly = deltaFPoly * deltaFPoly + deltaFSqPoly.changeBasis(("Array", self.basisM, self.basisN, self.basisN)) + + operator, _, _, collision = self.buildLinearEquations() + source = np.sum( + collision * deltaFSqPoly.coefficients[None, None, None, None, ...], + axis=(4, 5, 6, 7), + ) + + # Computing the correction from nonlinear terms + deltaNonlin = np.linalg.solve( + operator, np.reshape(source, source.size, order="C") + ) + deltaNonlinShape = ( + len(self.offEqParticles), + self.grid.M - 1, + self.grid.N - 1, + self.grid.N - 1, + ) + deltaNonlin = np.reshape(deltaNonlin, deltaNonlinShape, order="C") + deltaNonlinPoly = Polynomial( + coefficients=deltaNonlin, + grid=self.grid, + basis=("Array", self.basisM, self.basisN, self.basisN), + direction=("z", "z", "pz", "pp"), + endpoints=False, + ) + deltaNonlinPoly.changeBasis(("Array", "Cardinal", "Cardinal", "Cardinal")) + msqFull = np.array( [ particle.msqVacuum(self.background.fieldProfiles) for particle in particles ] ) - fieldPoly = Polynomial( - np.sum(self.background.fieldProfiles, axis=1), + + msqPoly = Polynomial( + msqFull, self.grid, - "Cardinal", + ("Array", "Cardinal"), "z", True, ) - dfielddChi = fieldPoly.derivative(0).coefficients[None, 1:-1, None, None] + dmsqdChi = msqPoly.derivative(axis=1).coefficients[:, 1:-1, None, None] # adding new axes, to make everything rank 3 like deltaF (z, pz, pp) # for fast multiplication of arrays, using numpy's broadcasting rules @@ -389,45 +605,24 @@ def checkLinearization( dpzdrz = dpzdrz[None, None, :, None] dppdrp = dppdrp[None, None, None, :] - # base integrand, for '00' - integrand = dfielddChi * dpzdrz * dppdrp * pp / (4 * np.pi**2 * energy) - - # The first criterion is to require that pressureOut/pressureEq is small - pressureOut = deltaFPoly.integrate((1, 2, 3), integrand).coefficients - pressureEq = fEqPoly.integrate((1, 2, 3), integrand).coefficients - deltaFCriterion = pressureOut / pressureEq - - # If criterion1 is large, we need C[deltaF]/L[deltaF] to be small - _, _, liouville, collision = self.buildLinearEquations() - collisionDeltaF = np.sum( - collision * deltaF[None, None, None, None, ...], axis=(4, 5, 6, 7) - ) - liouvilleDeltaF = np.sum( - liouville * deltaF[None, None, None, None, ...], axis=(4, 5, 6, 7) - ) - collisionDeltaFPoly = Polynomial( - collisionDeltaF, - self.grid, - ("Array", "Cardinal", "Cardinal", "Cardinal"), - ("z", "z", "pz", "pp"), - False, - ) - lioviilleDeltaFPoly = Polynomial( - liouvilleDeltaF, - self.grid, - ("Array", "Cardinal", "Cardinal", "Cardinal"), - ("z", "z", "pz", "pp"), - False, + dofs = np.array([particle.totalDOFs for particle in particles])[ + :, None, None, None + ] + integrand = dofs * dmsqdChi * dpzdrz * dppdrp * pp / (4 * np.pi**2 * energy) + + # Computing the pressure contributions of the equilibrium part, the linear + # out-of-equilibrium part and the first-order correction due to nonlinearities. + pressureEq = np.sum(fEqPoly.integrate((1, 2, 3), integrand).coefficients) + pressureDeltaF = np.sum(deltaFPoly.integrate((1, 2, 3), integrand).coefficients) + pressureNonlin = np.sum( + deltaNonlinPoly.integrate((1, 2, 3), integrand).coefficients ) - collisionDeltaFIntegrated = collisionDeltaFPoly.integrate( - (1, 2, 3), integrand - ).coefficients - liovilleDeltaFIntegrated = lioviilleDeltaFPoly.integrate( - (1, 2, 3), integrand - ).coefficients - collCriterion = collisionDeltaFIntegrated / liovilleDeltaFIntegrated - return deltaFCriterion, collCriterion + # Computing the 2 linearisation criteria + criterion1 = abs(pressureDeltaF / pressureEq) + criterion2 = abs(pressureNonlin / (pressureEq + pressureDeltaF)) + + return criterion1, criterion2 def buildLinearEquations( self, @@ -645,7 +840,7 @@ def loadCollisions(self, directoryPath: "pathlib.Path") -> None: self.basisN, self.offEqParticles, ) - logging.debug(f"Loaded collision data from directory {directoryPath}") + logging.debug("Loaded collision data from directory %s", directoryPath) except CollisionLoadError as e: raise diff --git a/src/WallGo/config.py b/src/WallGo/config.py index 9e5f8c1a..23a39dd2 100644 --- a/src/WallGo/config.py +++ b/src/WallGo/config.py @@ -119,6 +119,12 @@ class ConfigThermodynamics: uses the initial step size algorithm of :py:mod:`scipy.integrate.solve_ivp`. """ + interpolationDegree: int = 1 + """ + Degree of the splines used in FreeEnergy to interpolate the potential and its + derivatives. + """ + @dataclass class ConfigBoltzmannSolver: """ Holds the config of the BoltzmannSolver class. """ @@ -138,6 +144,12 @@ class ConfigBoltzmannSolver: WARNING: THIS CHANGES THE COLLISION TERMS WRT TO THEIR PHYSICAL VALUE. """ + truncationOption: str = 'AUTO' + """ Truncation option for spectral expansions. Can be 'NONE' for no + truncation, 'AUTO' to automatically detect if the spectral expansion + is converging and truncate if not, or 'THIRD' which always truncates + the last third. """ + @dataclass class Config: """ @@ -282,17 +294,24 @@ def loadConfigFromFile(self, filePath: str) -> None: self.configThermodynamics.phaseTracerFirstStep = None else: raise + if 'interpolationDegree' in keys: + self.configThermodynamics.interpolationDegree = parser.getint( + "Thermodynamics", + "interpolationDegree" + ) # Read the BoltzmannSolver configs if 'BoltzmannSolver' in parser.sections(): keys = parser['BoltzmannSolver'].keys() if 'collisionMultiplier' in keys: - self.configBoltzmannSolver.collisionMultiplier = parser.getfloat( - "BoltzmannSolver", - "collisionMultiplier") + self.configBoltzmannSolver.collisionMultiplier = parser.getfloat("BoltzmannSolver", "collisionMultiplier") if 'basisM' in keys: - self.configBoltzmannSolver.basisM = parser.get("BoltzmannSolver", - "basisM") + self.configBoltzmannSolver.basisM = parser.get( + "BoltzmannSolver", "basisM") if 'basisN' in keys: - self.configBoltzmannSolver.basisN = parser.get("BoltzmannSolver", - "basisN") + self.configBoltzmannSolver.basisN = parser.get( + "BoltzmannSolver", "basisN") + if 'truncationOption' in keys: + self.configBoltzmannSolver.truncationOption = parser.get( + "BoltzmannSolver", "truncationOption") + diff --git a/src/WallGo/containers.py b/src/WallGo/containers.py index 16b300d7..5289b8dd 100644 --- a/src/WallGo/containers.py +++ b/src/WallGo/containers.py @@ -4,6 +4,7 @@ from dataclasses import dataclass import numpy as np from .fields import Fields +from .freeEnergy import FreeEnergyValueType from .helpers import boostVelocity from .polynomial import Polynomial @@ -23,6 +24,56 @@ class PhaseInfo: """Temperature of transition.""" +@dataclass +class FreeEnergyArrays: + """Object containing temperatures, positions of the minimum and value of the effective potential. + """ + + temperatures: np.ndarray[float] # 1D array + """Array of temperatures.""" + + freeEnergyList: np.ndarray[FreeEnergyValueType] # 1D array of FreeEnergyValueType objects + """Array of field(s) and potential value at the minimum.""" + + allowedDiscrepancy: float | None + """Allowed discrepancy between the effective potential at the minimum and the user-provided value""" + + def __init__( + self, + temperatures: np.ndarray[float], + minimumList: np.ndarray[float], + potentialEffList: np.ndarray[float], + allowedDiscrepancy: float | None = None, + ) -> None: + """Initialisation of FreeEnergyArrays, based on passing 3 arrays + and a float. """ + temperatures = np.asarray(temperatures) + minimumList = np.asarray(minimumList) + potentialEffList = np.asarray(potentialEffList) + + if temperatures.ndim != 1: + raise ValueError("temperatures must be a 1D array.") + if temperatures.shape[0] != minimumList.shape[0]: + raise ValueError( + "The temperatures and minimumList must have the same length." + ) + if temperatures.shape[0] != potentialEffList.shape[0]: + raise ValueError( + "The temperatures and potentialEffList must have the same length." + ) + if allowedDiscrepancy is not None and allowedDiscrepancy < 0: + raise ValueError("allowedDiscrepancy must not be negative.") + + self.temperatures = temperatures + self.freeEnergyList = FreeEnergyValueType.fromArray( + np.concatenate( + (minimumList, potentialEffList[:, np.newaxis]), axis=1, + ) + ) + self.allowedDiscrepancy = allowedDiscrepancy + + + @dataclass class WallParams: """ diff --git a/src/WallGo/effectivePotential.py b/src/WallGo/effectivePotential.py index 838af60d..b142a212 100644 --- a/src/WallGo/effectivePotential.py +++ b/src/WallGo/effectivePotential.py @@ -120,7 +120,8 @@ def getInherentRelativeError(self) -> float: return self.effectivePotentialError - def findLocalMinimum(self, initialGuess: Fields, temperature: npt.ArrayLike, tol: float = None) -> Tuple[Fields, np.ndarray]: + def findLocalMinimum(self, initialGuess: Fields, temperature: npt.ArrayLike, + tol: float = None, method: str|None = None) -> Tuple[Fields, np.ndarray]: """ Finds a local minimum starting from a given initial configuration of background fields. Feel free to override this if your model requires more delicate minimization. @@ -159,7 +160,7 @@ def evaluateWrapper(fieldArray: np.ndarray): guess = guesses.getFieldPoint(i) - res = scipy.optimize.minimize(evaluateWrapper, guess, tol=tol) + res = scipy.optimize.minimize(evaluateWrapper, guess, tol=tol, method=method) resLocation[i] = res.x resValue[i] = res.fun diff --git a/src/WallGo/equationOfMotion.py b/src/WallGo/equationOfMotion.py index ffc9bd23..bc669eda 100644 --- a/src/WallGo/equationOfMotion.py +++ b/src/WallGo/equationOfMotion.py @@ -53,6 +53,7 @@ def __init__( errTol: float = 1e-3, maxIterations: int = 10, pressRelErrTol: float = 0.3679, + # pylint: disable=too-many-arguments ): """ Initialization @@ -136,6 +137,11 @@ def __init__( ## Flag to detect if we were able to find the pressure self.successWallPressure = True + ## Setup lists used to estimate the pressure derivative + self.listVelocity = [] + self.listPressure = [] + self.listPressureError = [] + def findWallVelocityDeflagrationHybrid( self, wallThicknessIni: float | None = None ) -> WallGoResults: @@ -176,7 +182,7 @@ def findWallVelocityDeflagrationHybrid( vmax = min(self.hydrodynamics.vJ, self.hydrodynamics.fastestDeflag()) if vmax < self.hydrodynamics.vJ and ( - self.hydrodynamics.doesPhaseTraceLimitvmax[0] + self.hydrodynamics.doesPhaseTraceLimitvmax[0] or self.hydrodynamics.doesPhaseTraceLimitvmax[1] ): logging.warning( @@ -428,6 +434,11 @@ def solveWall( Data class containing results. """ + ## Reset lists used to estimate the pressure derivative + self.listVelocity = [] + self.listPressure = [] + self.listPressureError = [] + results = WallGoResults() results.hasOutOfEquilibrium = self.includeOffEq @@ -441,6 +452,8 @@ def solveWall( boltzmannResultsMax, boltzmannBackgroundMax, hydroResultsMax, + emViolationT30Max, + emViolationT33Max, ) = self.wallPressure(wallVelocityMax, wallParamsGuess) else: ( @@ -449,6 +462,8 @@ def solveWall( boltzmannResultsMax, boltzmannBackgroundMax, hydroResultsMax, + emViolationT30Max, + emViolationT33Max, ) = wallPressureResultsMax # also getting the LTE results @@ -458,12 +473,13 @@ def solveWall( # hybrid solution if pressureMax < 0: logging.info("Maximum pressure on wall is negative!") - logging.info(f"{pressureMax=} {wallParamsMax=}") + logging.info("pressureMax=%s wallParamsMax=%s", pressureMax, wallParamsMax) results.setWallVelocities(None, None, wallVelocityLTE) results.setWallParams(wallParamsMax) results.setHydroResults(hydroResultsMax) results.setBoltzmannBackground(boltzmannBackgroundMax) results.setBoltzmannResults(boltzmannResultsMax) + results.setViolationOfEMConservation((emViolationT30Max, emViolationT33Max)) results.setSuccessState( True, ESolutionType.RUNAWAY, @@ -480,6 +496,8 @@ def solveWall( boltzmannResultsMin, boltzmannBackgroundMin, hydroResultsMin, + emViolationT30Min, + emViolationT33Min, ) = self.wallPressure(wallVelocityMin, wallParamsGuess) else: ( @@ -488,6 +506,8 @@ def solveWall( boltzmannResultsMin, boltzmannBackgroundMin, hydroResultsMin, + emViolationT30Min, + emViolationT33Min, ) = wallPressureResultsMin while pressureMin > 0: @@ -505,6 +525,9 @@ def solveWall( results.setHydroResults(hydroResultsMin) results.setBoltzmannBackground(boltzmannBackgroundMin) results.setBoltzmannResults(boltzmannResultsMin) + results.setViolationOfEMConservation( + (emViolationT30Min, emViolationT33Min) + ) results.setSuccessState( False, ESolutionType.ERROR, @@ -518,6 +541,8 @@ def solveWall( boltzmannResultsMin, boltzmannBackgroundMin, hydroResultsMin, + emViolationT30Min, + emViolationT33Min, ) = self.wallPressure(wallVelocityMin, wallParamsGuess) self.pressAbsErrTol = ( @@ -576,31 +601,73 @@ def pressureWrapper(vw: float) -> float: # pylint: disable=invalid-name boltzmannResults, boltzmannBackground, hydroResults, + emViolationT30, + emViolationT33, ) = self.wallPressure( wallVelocity, newWallParams, boltzmannResultsInput=newBoltzmannResults ) + eomResidual = self.estimateTanhError( + wallParams, boltzmannResults, boltzmannBackground, hydroResults + ) + # minimum possible error in the wall speed - wallVelocityMinError = self.errTol * optimizeResult.root + wallVelocityMinError = self.errTol - # estimating errors from truncation and comparison to finite differences if self.includeOffEq: + # Computing the linearisation criteria + criterion1, criterion2 = self.boltzmannSolver.checkLinearization( + boltzmannResults.deltaF + ) + boltzmannResults.linearizationCriterion1 = criterion1 + boltzmannResults.linearizationCriterion2 = criterion2 + + # Computing the out-of-equilibrium pressure to get the absolute error + vevLowT = boltzmannBackground.fieldProfiles.getFieldPoint(0) + vevHighT = boltzmannBackground.fieldProfiles.getFieldPoint(-1) + fields, dPhidz = self.wallProfile( + self.grid.xiValues, vevLowT, vevHighT, wallParams + ) + dVout = ( + np.sum( + [ + particle.totalDOFs + * particle.msqDerivative(fields) + * boltzmannResults.Deltas.Delta00.coefficients[i, :, None] + for i, particle in enumerate(self.particles) + ], + axis=0, + ) + / 2 + ) + + dVoutdz = np.sum(np.array(dVout * dPhidz), axis=1) + + # Create a Polynomial object to represent dVdz. Will be used to integrate. + dVoutdzPoly = Polynomial(dVoutdz, self.grid) + + dzdchi, _, _ = self.grid.getCompactificationDerivatives() + offEquilPressureScale = np.abs(dVoutdzPoly.integrate(weight=-dzdchi)) + + # Compute the pressure derivative + pressureDerivative = self.estimatePressureDerivative(wallVelocity) + + # estimating errors from truncation and comparison to finite differences finiteDifferenceBoltzmannResults = self.getBoltzmannFiniteDifference() - # assuming nonequilibrium errors proportional to deviation from LTE - wallVelocityDeltaLTE = abs(wallVelocity - wallVelocityLTE) # the truncation error in the spectral method within Boltzmann - wallVelocityTruncationError = ( - boltzmannResults.truncationError * wallVelocityDeltaLTE + wallVelocityTruncationError = abs( + boltzmannResults.truncationError + * offEquilPressureScale + / pressureDerivative ) # the deviation from the finite difference method within Boltzmann delta00 = boltzmannResults.Deltas.Delta00.coefficients[0] delta00FD = finiteDifferenceBoltzmannResults.Deltas.Delta00.coefficients[0] errorFD = np.linalg.norm(delta00 - delta00FD) / np.linalg.norm(delta00) - wallVelocityDerivativeError = errorFD * wallVelocityDeltaLTE # if truncation waringin large, raise a warning if ( - wallVelocityTruncationError > wallVelocityDerivativeError + boltzmannResults.truncationError > errorFD and wallVelocityTruncationError > self.errTol ): warnings.warn("Truncation error large, increase N or M", RuntimeWarning) @@ -626,6 +693,8 @@ def pressureWrapper(vw: float) -> float: # pylint: disable=invalid-name results.setBoltzmannBackground(boltzmannBackground) results.setBoltzmannResults(boltzmannResults) results.setFiniteDifferenceBoltzmannResults(finiteDifferenceBoltzmannResults) + results.setViolationOfEMConservation((emViolationT30, emViolationT33)) + results.eomResidual = eomResidual # Set the message if not self.successTemperatureProfile: @@ -695,7 +764,15 @@ def wallPressure( atol: float | None = None, rtol: float | None = None, boltzmannResultsInput: BoltzmannResults | None = None, - ) -> tuple[float, WallParams, BoltzmannResults, BoltzmannBackground, HydroResults]: + ) -> tuple[ + float, + WallParams, + BoltzmannResults, + BoltzmannBackground, + HydroResults, + float, + float, + ]: """ Computes the total pressure on the wall by finding the tanh profile that minimizes the action. Can use two different iteration algorithms @@ -737,7 +814,12 @@ def wallPressure( hydroResults : HydroResults HydroResults object containing the solution obtained from Hydrodynamics. Only returned if returnExtras is True - + emViolationAfter[0] : float + Violation of energy-momentum conservation in T30 after solving the + Boltzmann equation + emViolationAfter[1] : float + Violation of energy-momentum conservation in T33 after solving the + Boltzmann equation """ if atol is None: @@ -747,11 +829,11 @@ def wallPressure( self.successWallPressure = True - improveConvergence = self.forceImproveConvergence + cautious = self.forceImproveConvergence if wallVelocity > self.hydrodynamics.vJ: - improveConvergence = True + cautious = True - logging.info(f"------------- Trying {wallVelocity=:g} -------------") + logging.info("------------- Trying wallVelocity=%g -------------", wallVelocity) # Initialize the different data class objects and arrays zeroPoly = Polynomial( @@ -781,8 +863,8 @@ def wallPressure( deltaF=deltaF, Deltas=offEquilDeltas, truncationError=0.0, - linearizationCriterion1=np.zeros(len(self.particles)), - linearizationCriterion2=np.zeros(len(self.particles)), + truncatedTail=(False, False, False), + spectralPeaks=(0, 0, 0), ) else: boltzmannResults = boltzmannResultsInput @@ -823,6 +905,8 @@ def wallPressure( wallParams, boltzmannResults, boltzmannBackground, + emViolationBefore, + emViolationAfter, ) = self._intermediatePressureResults( wallParams, vevLowT, @@ -857,11 +941,34 @@ def wallPressure( multiplier = 1.0 i = 0 - logging.debug( - f"{'pressure':>12s} {'error':>12s} {'errorSolver':>12s} {'errTol':>12s} {'cautious':>12s} {'multiplier':>12s}" - ) + if self.includeOffEq: + logging.debug( + "%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s", + "pressure", + "error", + "errorSolver", + "errTol", + "cautious", + "multiplier", + "dT30Before", + "dT30After", + "spectralPeak", + "truncatedTail", + ) + else: + logging.debug( + "%12s %12s %12s %12s %12s %12s %12s", + "pressure", + "error", + "errorSolver", + "errTol", + "cautious", + "multiplier", + "dT30Before", + ) + while True: - if improveConvergence: + if cautious: # Use the improved algorithm (which converges better but slowly) ( pressure, @@ -869,6 +976,8 @@ def wallPressure( boltzmannResults, boltzmannBackground, errorSolver, + emViolationBefore, + emViolationAfter, ) = self._getNextPressure( pressure, wallParams, @@ -890,6 +999,8 @@ def wallPressure( wallParams, boltzmannResults, boltzmannBackground, + emViolationBefore, + emViolationAfter, ) = self._intermediatePressureResults( wallParams, vevLowT, @@ -910,12 +1021,34 @@ def wallPressure( error = np.abs(pressures[-1] - pressures[-2]) errTol = np.maximum(rtol * np.abs(pressure), atol) * multiplier - logging.debug( - f"{pressure:>12g} {error:>12g} {errorSolver:>12g} {errTol:>12g} {improveConvergence:>12} {multiplier:>12g}" - ) + if self.includeOffEq: + logging.debug( + "%12g %12g %12g %12g %12r %12g %12g %12g %12r %12r", + pressure, + error, + errorSolver, + errTol, + int(cautious), + multiplier, + emViolationBefore[0], + emViolationAfter[0], + tuple(int(s) for s in boltzmannResults.spectralPeaks), + tuple(int(t) for t in boltzmannResults.truncatedTail), + ) + else: + logging.debug( + "%12g %12g %12g %12g %12r %12g %12g", + pressure, + error, + errorSolver, + errTol, + int(cautious), + multiplier, + emViolationBefore[0], + ) i += 1 - if error < errTol or (errorSolver < errTol and improveConvergence): + if error < errTol or (errorSolver < errTol and cautious): """ Even if two consecutive call to _getNextPressure() give similar pressures, it is possible that the internal calls made to @@ -952,17 +1085,23 @@ def wallPressure( if len(pressures) > 2: if error > abs(pressures[-2] - pressures[-3]) / 1.5: # If the error decreases too slowly, use the improved algorithm - improveConvergence = True + cautious = True logging.info(f"Final {pressure=:g}") logging.debug(f"Final {wallParams=}") - + + self.listVelocity.append(wallVelocity) + self.listPressure.append(pressure) + self.listPressureError.append(max(error, rtol * np.abs(pressure), atol)) + return ( pressure, wallParams, boltzmannResults, boltzmannBackground, hydroResults, + emViolationAfter[0], + emViolationAfter[1], ) def _getNextPressure( @@ -987,13 +1126,20 @@ def _getNextPressure( Boltzmann results each time. If the iterations overshot the true solution (the two pressure updates go in opposite directions), uses linear interpolation to find a better estimate of the true solution. This function is - called only when improveConvergence=True in wallPressure(). + called only when cautious=True in wallPressure(). """ + + # The different pressure_i, wallParams_i and boltzmannResults_i correspond to + # different iterations that the solver use to determine if it overshoots the + # true solution. If it does, it uses linear interpolation to find a better + # solution. ( pressure2, wallParams2, boltzmannResults2, _, + _, + _, ) = self._intermediatePressureResults( wallParams1, vevLowT, @@ -1013,6 +1159,8 @@ def _getNextPressure( wallParams3, boltzmannResults3, boltzmannBackground3, + emViolationBefore, + emViolationAfter, ) = self._intermediatePressureResults( wallParams2, vevLowT, @@ -1032,7 +1180,15 @@ def _getNextPressure( ## last update go in the same direction), returns the last iteration. if (pressure3 - pressure2) * (pressure2 - pressure1) >= 0: err = abs(pressure3 - pressure2) - return pressure3, wallParams3, boltzmannResults3, boltzmannBackground3, err + return ( + pressure3, + wallParams3, + boltzmannResults3, + boltzmannBackground3, + err, + emViolationBefore, + emViolationAfter, + ) ## If the last iteration overshot, uses linear interpolation to find a ## better estimate of the true solution. @@ -1042,6 +1198,8 @@ def _getNextPressure( wallParams4, boltzmannResults4, boltzmannBackground4, + emViolationBefore, + emViolationAfter, ) = self._intermediatePressureResults( wallParams1 + (wallParams2 - wallParams1) * interpPoint, vevLowT, @@ -1057,7 +1215,15 @@ def _getNextPressure( multiplier, ) err = abs(pressure4 - pressure2) - return pressure4, wallParams4, boltzmannResults4, boltzmannBackground4, err + return ( + pressure4, + wallParams4, + boltzmannResults4, + boltzmannBackground4, + err, + emViolationBefore, + emViolationAfter, + ) def _intermediatePressureResults( self, @@ -1073,7 +1239,14 @@ def _intermediatePressureResults( temperatureProfileInput: np.ndarray | None = None, velocityProfileInput: np.ndarray | None = None, multiplier: float = 1.0, - ) -> tuple[float, WallParams, BoltzmannResults, BoltzmannBackground]: + ) -> tuple[ + float, + WallParams, + BoltzmannResults, + BoltzmannBackground, + tuple[float, float], + tuple[float, float], + ]: """ Performs one step of the iteration procedure to update the pressure, wall parameters and Boltzmann solution. This is done by first solving @@ -1113,6 +1286,20 @@ def _intermediatePressureResults( temperatureProfile = temperatureProfileInput velocityProfile = velocityProfileInput + # Compute the violation of energy-momentum conservation before solving + # the Boltzmann equation + violationOfEMConservationBefore = self.violationOfEMConservation( + c1, + c2, + velocityMid, + fields, + dfieldsdz, + boltzmannResults.Deltas, + temperatureProfile, + velocityProfile, + max(wallParams.widths), + ) + ## Prepare a new background for Boltzmann TWithEndpoints: np.ndarray = np.concatenate( (np.array([Tminus]), np.array(temperatureProfile), np.array([Tplus])) @@ -1187,6 +1374,20 @@ def actionWrapper( ) dVdPhi = self.thermo.effectivePotential.derivField(fields, temperatureProfile) + # Compute the violation of energy-momentum conservation after + # solving the Boltzmann equation + violationOfEMConservationAfter = self.violationOfEMConservation( + c1, + c2, + velocityMid, + fields, + dPhidz, + boltzmannResults.Deltas, + temperatureProfile, + velocityProfile, + max(wallParams.widths), + ) + # Out-of-equilibrium term of the EOM dVout = ( np.sum( @@ -1212,7 +1413,14 @@ def actionWrapper( dzdchi, _, _ = self.grid.getCompactificationDerivatives() pressure = eomPoly.integrate(weight=-dzdchi) - return pressure, wallParams, boltzmannResults, boltzmannBackground + return ( + pressure, + wallParams, + boltzmannResults, + boltzmannBackground, + violationOfEMConservationBefore, + violationOfEMConservationAfter, + ) def _toWallParams(self, wallArray: np.ndarray) -> WallParams: offsets: np.ndarray = np.concatenate( @@ -1344,6 +1552,133 @@ def action( return float(U + K) + def estimateTanhError( + self, + wallParams: WallParams, + boltzmannResults: BoltzmannResults, + boltzmannBackground: BoltzmannBackground, + hydroResults: HydroResults, + ) -> np.ndarray: + r""" + Estimates the EOM error due to the tanh ansatz. It is estimated by the integral + + .. math:: \sqrt{\Delta[\mathrm{EOM}^2]/|\mathrm{EOM}^2|}, + + with + + .. math:: \\Delta[\\mathrm{EOM}^2]=\\int\\! dz\\, (-\\partial_z^2 \\phi+ + \\partial V_{\\mathrm{eq}}/ \\partial \\phi+ \\partial V_{\\mathrm{out}}/ \\partial \\phi )^2 + + and + + .. math:: |\\mathrm{EOM}^2|=\\int\\! dz\\, [(\\partial_z^2 \\phi)^2+ + (\\partial V_{\\mathrm{eq}}/ \\partial \\phi)^2+ (\\partial V_{\\mathrm{out}}/ \\partial \\phi)^2]. + + """ + Tminus = hydroResults.temperatureMinus + Tplus = hydroResults.temperaturePlus + + # Positions of the phases + TminusEval = max( + min(Tminus, self.thermo.freeEnergyLow.interpolationRangeMax()), + self.thermo.freeEnergyLow.interpolationRangeMin(), + ) + TplusEval = max( + min(Tplus, self.thermo.freeEnergyHigh.interpolationRangeMax()), + self.thermo.freeEnergyHigh.interpolationRangeMin(), + ) + vevLowT = self.thermo.freeEnergyLow(TminusEval).fieldsAtMinimum + vevHighT = self.thermo.freeEnergyHigh(TplusEval).fieldsAtMinimum + + temperatureProfile = boltzmannBackground.temperatureProfile[1:-1] + + z = self.grid.xiValues + fields = self.wallProfile(z, vevLowT, vevHighT, wallParams)[0] + with warnings.catch_warnings(): + # overflow here is benign, as just gives zero + warnings.filterwarnings("ignore", message="overflow encountered in *") + d2FieldsDz2 = -( + (vevHighT - vevLowT) + * np.tanh(z[:, None] / wallParams.widths[None, :] + wallParams.offsets) + / np.cosh(z[:, None] / wallParams.widths[None, :] + wallParams.offsets) ** 2 + / wallParams.widths**2 + ) + + dVdPhi = self.thermo.effectivePotential.derivField(fields, temperatureProfile) + + # Out-of-equilibrium term of the EOM + dVout = ( + np.sum( + [ + particle.totalDOFs + * particle.msqDerivative(fields) + * boltzmannResults.Deltas.Delta00.coefficients[i, :, None] + for i, particle in enumerate(self.particles) + ], + axis=0, + ) + / 2 + ) + + eomSq = (-d2FieldsDz2 + dVdPhi + dVout) ** 2 + eomSqScale = d2FieldsDz2**2 + dVdPhi**2 + dVout**2 + + eomSqPoly = Polynomial(eomSq, self.grid, basis=("Cardinal", "Array")) + eomSqScalePoly = Polynomial(eomSqScale, self.grid, basis=("Cardinal", "Array")) + dzdchi, _, _ = self.grid.getCompactificationDerivatives() + eomSqResidual = eomSqPoly.integrate(axis=0, weight=dzdchi[:, None]) + eomSqScaleIntegrated = eomSqScalePoly.integrate(axis=0, weight=dzdchi[:, None]) + + return eomSqResidual.coefficients / eomSqScaleIntegrated.coefficients + + def estimatePressureDerivative(self, wallVelocity: float) -> float: + """ + Estimates the derivative of the preessure with respect to the wall velocity from + a least square fit of the computed pressure to a line. Must have run + wallPressure at velocities close to wallVelocity before calling this function. + + Parameters + ---------- + wallVelocity : float + Wall velocity. + + Returns + ------- + float + Derivative of the pressure at wallVelocity. + + """ + # Number of pressure points + nbrPressure = len(self.listPressure) + + assert ( + len(self.listPressureError) == len(self.listVelocity) == nbrPressure >= 2 + ), """The lists listVelocity, listPressure, + listPressureError must have the same length and + contain at least two elements.""" + + velocityErrorScale = self.errTol * wallVelocity + pressures = np.array(self.listPressure) + velocityDiff = np.array(self.listVelocity) - wallVelocity + # Farter points are exponentially suppressed to make sure they don't impact the + # estimate too much. + weightMatrix = np.diag( + np.exp(-np.abs(velocityDiff / velocityErrorScale)) + / np.array(self.listPressureError) ** 2 + ) + aMatrix = np.ones((nbrPressure, 2)) + aMatrix[:, 1] = velocityDiff + + # Computes the derivative by fitting the pressure to a line + derivative = ( + np.linalg.inv(aMatrix.T @ weightMatrix @ aMatrix) + @ aMatrix.T + @ weightMatrix + @ pressures + )[1] + + return derivative + def wallProfile( self, z: np.ndarray, # pylint: disable=invalid-name @@ -1382,7 +1717,8 @@ def wallProfile( zL = z[:, None] / wallParams.widths[None, :] # pylint: disable=invalid-name with warnings.catch_warnings(): - warnings.simplefilter("ignore") + # overflow here is benign, as just gives zero + warnings.filterwarnings("ignore", message="overflow encountered in *") fields = vevLowT + 0.5 * (vevHighT - vevLowT) * ( 1 + np.tanh(zL + wallParams.offsets) ) @@ -1585,7 +1921,7 @@ def plasmaVelocity( Plasma velocity. """ - # Need enthalpy ouside a free-energy minimum (eq .(12) in arXiv:2204.13120v1) + # Need enthalpy outside a free-energy minimum (eq .(12) in arXiv:2204.13120v1) enthalpy = -T * self.thermo.effectivePotential.derivT(fields, T) return float((-enthalpy + np.sqrt(4 * s1**2 + enthalpy**2)) / (2 * s1)) @@ -1620,7 +1956,7 @@ def temperatureProfileEqLHS( LHS of Eq. (20) of arXiv:2204.13120v1. """ - # Need enthalpy ouside a free-energy minimum (eq (12) in the ref.) + # Need enthalpy outside a free-energy minimum (eq (12) in the ref.) enthalpy = -T * self.thermo.effectivePotential.derivT(fields, T) kineticTerm = 0.5 * np.sum(dPhidz**2).view(np.ndarray) @@ -1766,3 +2102,153 @@ def getBoltzmannFiniteDifference(self) -> BoltzmannResults: boltzmannSolverFiniteDifference.collisionArray.changeBasis("Cardinal") # now computing results return boltzmannSolverFiniteDifference.getDeltas() + + def violationOfEMConservation( + self, + c1: float, # pylint: disable=invalid-name + c2: float, # pylint: disable=invalid-name + velocityMid: float, + fields: Fields, + dPhidz: Fields, + offEquilDeltas: BoltzmannDeltas, + temperatureProfile: np.ndarray, + velocityProfile: np.ndarray, + wallThickness: float, + ) -> Tuple[float, float]: + r""" + Determines the RMS (along the grid) of the residual of the + energy-momentum equations (18) of arXiv:2204.13120v1. + + Parameters + ---------- + index : int + Index of the grid point. + c1 : float + Value of the :math:`T^{30}` component of the energy-momentum tensor. + c2 : float + Value of the :math:`T^{33}` component of the energy-momentum tensor. + velocityMid : float + Midpoint of plasma velocity in the wall frame, :math:`(v_+ + v_-)/2`. + fields : FieldPoint + Scalar field profile. + dPhidz : FieldPoint + Derivative with respect to the position of the scalar field profile. + offEquilDeltas : BoltzmannDeltas + BoltzmannDeltas object containing the off-equilibrium Delta functions + temperatureProfile: np.ndarray + Plasma temperature profile at the grid points. + velocityProfile: np.ndarray + Plasma velocity profile at the grid points. + wallThickness: float + Thickness of the wall, used to normalize the violations. + + Returns + ------- + violationEM30, violationEM33 : (float, float) + Violation of energy-momentum conservation in T03 and T33 integrated over + the grid, normalized by the wall thickness. + + """ + + violationT30sq = np.zeros(len(self.grid.xiValues)) + violationT33sq = np.zeros(len(self.grid.xiValues)) + + for index in range(len(self.grid.xiValues)): + vt30, vt33 = self.violationEMPoint( + index, + c1, + c2, + velocityMid, + fields.getFieldPoint(index), + dPhidz.getFieldPoint(index), + offEquilDeltas, + temperatureProfile[index], + velocityProfile[index], + ) + + violationT30sq[index] = ( + (np.asarray(vt30).item()) ** 2 + if isinstance(vt30, np.ndarray) + else vt30**2 + ) + violationT33sq[index] = ( + (np.asarray(vt33).item()) ** 2 + if isinstance(vt33, np.ndarray) + else vt33**2 + ) + + T30Poly = Polynomial(violationT30sq, self.grid) + T33Poly = Polynomial(violationT33sq, self.grid) + dzdchi, _, _ = self.grid.getCompactificationDerivatives() + violationT30sqIntegrated = T30Poly.integrate(weight=dzdchi) + violationT33sqIntegrated = T33Poly.integrate(weight=dzdchi) + + return ( + np.sqrt(violationT30sqIntegrated) / wallThickness, + np.sqrt(violationT33sqIntegrated) / wallThickness, + ) + + def violationEMPoint( + self, + index: int, + c1: float, # pylint: disable=invalid-name + c2: float, # pylint: disable=invalid-name + velocityMid: float, + fields: FieldPoint, + dPhidz: FieldPoint, + offEquilDeltas: BoltzmannDeltas, + T: float, + v: float, # pylint: disable=invalid-name + ) -> Tuple[float, float]: + r""" + Determines the residual of the energy-momentum equations (18) of + arXiv:2204.13120v1 locally. + + Parameters + ---------- + index : int + Index of the grid point. + c1 : float + Value of the :math:`T^{30}` component of the energy-momentum tensor. + c2 : float + Value of the :math:`T^{33}` component of the energy-momentum tensor. + velocityMid : float + Midpoint of plasma velocity in the wall frame, :math:`(v_+ + v_-)/2`. + fields : FieldPoint + Scalar field profile. + dPhidz : FieldPoint + Derivative with respect to the position of the scalar field profile. + offEquilDeltas : BoltzmannDeltas + BoltzmannDeltas object containing the off-equilibrium Delta functions + T: float + Plasma temperature at the point grid.xiValues[index]. + v: float + Plasma velocity at the point grid.xiValues[index]. + + Returns + ------- + violationEM30, violationEM33 : float + Violation of energy-momentum conservation in T03 and T33 at + the point grid.xiValues[index]. + + """ + + # Computing the out-of-equilibrium part of the energy-momentum tensor + Tout30, Tout33 = self.deltaToTmunu(index, fields, velocityMid, offEquilDeltas) + + # Need enthalpy ouside a free-energy minimum (eq .(12) in arXiv:2204.13120v1) + enthalpy = -T * self.thermo.effectivePotential.derivT(fields, T) + + # Kinetic term + kineticTerm = 0.5 * np.sum(dPhidz**2).view(np.ndarray) + + ## eff potential at this field point and temperature. NEEDS the T-dep constant + veff = self.thermo.effectivePotential.evaluate(fields, T) + + violationEM30 = (enthalpy * v / (1 - v**2) + Tout30 - c1) / c1 + + violationEM33 = ( + kineticTerm - veff + enthalpy * v**2 / (1 - v**2) + Tout33 - c2 + ) / c2 + + return (violationEM30, violationEM33) diff --git a/src/WallGo/exceptions.py b/src/WallGo/exceptions.py index 4c9fcae1..e8295d60 100644 --- a/src/WallGo/exceptions.py +++ b/src/WallGo/exceptions.py @@ -3,7 +3,6 @@ """ import typing -from .containers import PhaseInfo class WallGoError(Exception): diff --git a/src/WallGo/freeEnergy.py b/src/WallGo/freeEnergy.py index 7cce8f8c..f525e403 100644 --- a/src/WallGo/freeEnergy.py +++ b/src/WallGo/freeEnergy.py @@ -2,21 +2,25 @@ 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 import scipy.integrate as scipyint import scipy.linalg as scipylinalg +from typing import Union + +# from .containers import FreeEnergyArrays +from .effectivePotential import EffectivePotential +from .exceptions import WallGoError +from .fields import FieldPoint, Fields from .interpolatableFunction import ( InterpolatableFunction, EExtrapolationType, inputType, outputType, ) -from .effectivePotential import EffectivePotential -from .exceptions import WallGoError -from .fields import FieldPoint, Fields @dataclass @@ -163,9 +167,9 @@ def __call__( "Trying to evaluate FreeEnergy outside of its range of existence" ) raise WallGoError( - """\n Trying to evaluate FreeEnergy outside of its allowed range, + """\n Trying to evaluate FreeEnergy outside of its allowed range, try increasing/decreasing Tmax/Tmin.""" - ) + ) def _functionImplementation(self, temperature: inputType | float) -> outputType: """ @@ -246,6 +250,7 @@ def tracePhase( 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 + interpolationDegree: int = 1, ) -> None: r"""Traces minimum of potential @@ -254,7 +259,7 @@ def tracePhase( .. math:: \frac{\partial^2 V^\text{eff}}{\partial \phi_i \partial \phi_j}\bigg|_{\phi=\phi^\text{min}} \frac{\partial \phi^\text{min}_j}{\partial T} + \frac{\partial^2 V^\text{eff}}{\partial \phi_i \partial T}\bigg|_{\phi=\phi^\text{min}} = 0, - + starting from a solution at the starting temperature. It uses `scipy.integrate.solve_ivp` to solve the problem. Stops if a mass squared goes through zero. Parameters @@ -273,6 +278,9 @@ def tracePhase( 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` + interpolationDegree : int, optional + Degree of the splines used in FreeEnergy to interpolate the potential and + its derivatives. Default is 1. """ # make sure the initial conditions are extra accurate extraTol = 0.01 * rTol @@ -296,17 +304,14 @@ def odeFunction(temperature: float, field: np.ndarray) -> np.ndarray: ) return np.asarray(scipylinalg.solve(hess, -dgraddT, assume_a="sym")) - # finding some sensible mass scales - ddVT0 = self.effectivePotential.deriv2Field2(phase0, T0) - eigsT0 = np.linalg.eigvalsh(ddVT0) - # mass_scale_T0 = np.mean(eigs_T0) - # min_mass_scale = rTol * mass_scale_T0 - # mass_hierarchy_T0 = min(eigs_T0) / max(eigs_T0) - # min_hierarchy = rTol * mass_hierarchy_T0 + # compute all the second derivatives at the beginning of phase tracing + d2Vdphi2, d2VdphidT, d2VdT2 = self.effectivePotential.allSecondDerivatives( + phase0, T0) + eigsT0 = np.linalg.eigvalsh(d2Vdphi2) # checking stable phase at initial temperature assert ( - min(eigsT0) * max(eigsT0) > 0 + min(eigsT0) > 0 ), "tracePhase error: unstable at starting temperature" def spinodalEvent(temperature: float, field: np.ndarray) -> float: @@ -316,12 +321,17 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: d2V = self.effectivePotential.deriv2Field2(FieldPoint(field), temperature) eigs = scipylinalg.eigvalsh(d2V) return float(min(eigs)) + # arrays to store results TList = np.full(1, T0) fieldList = np.full((1, phase0.numFields()), Fields((phase0,))) potentialEffList = np.full((1, 1), [potential0]) - + dVdTList = np.full((1,1), [self.effectivePotential.derivT(phase0, T0)]) + dphidT = -np.linalg.inv(d2Vdphi2)@d2VdphidT + d2VdT2List = np.full((1,1), [d2VdT2+dphidT@d2VdphidT]) + dphidTList = np.full((1, phase0.numFields()), + Fields((dphidT,))) # maximum temperature range TMin = max(self.minPossibleTemperature[0], TMin) TMax = min(self.maxPossibleTemperature[0], TMax) @@ -348,6 +358,9 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: while ode.status == "running": try: ode.step() + # check if all the eigenvalues of the hessian are positive + if spinodalEvent(ode.t, ode.y) <= 0: + break except RuntimeWarning as error: logging.error(error.args[0] + f" at T={ode.t}") break @@ -358,9 +371,7 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: tol=rTol, ) ode.y = phaset[0] - if spinodalEvent(ode.t, ode.y) <= 0: - break - if not paranoid: + else: # check if extremum is still accurate dVt = self.effectivePotential.derivField(Fields((ode.y)), ode.t) err = np.linalg.norm(dVt) / T0**3 @@ -378,28 +389,70 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: potentialEffT = np.asarray( self.effectivePotential.evaluate(Fields((ode.y)), ode.t) ) + + # Computing all the derivatives along the whole phase tracing + dVdT = self.effectivePotential.derivT(Fields((ode.y)), ode.t) + (d2Vdphi2, + d2VdphidT, + d2VdT2) = self.effectivePotential.allSecondDerivatives( + FieldPoint(ode.y), ode.t) + # check if step size is still okay to continue if ode.step_size < 1e-16 * T0 or ( TList.size > 0 and ode.t == TList[-1] ): logging.warning( - f"Step size {ode.step_size} shrunk too small at T={ode.t}, " - f"vev={ode.y}" + "Step size %g shrunk too small at T=%g, vev=%g", + ode.step_size, + ode.t, + ode.y, ) break + + dphidT = -np.linalg.inv(d2Vdphi2)@d2VdphidT + D2VDT2 = d2VdT2+dphidT@d2VdphidT + + # check that sound speed square is still positive + csq = dVdT/D2VDT2/ode.t + if csq < 0: + break + # Check if 2 methods for computing the 2nd derivative disagree by more + # than a factor of 2. This would indicate a discontinuity caused by a + # phase disappearing. + if TList.size >= 2: + # The first method uses the last value stored in d2VdT2List, which + # computes the total derivative in terms of the partial derivatives + # of V and the field phi. This is the more accurate method. + # The second method takes the finite derivative of dVdT, which + # should break down when the phase disappears because there is a + # discontinuity. + if d2VdT2List[-1,0]*(ode.t-TList[-2])/(dVdT-dVdTList[-2,0]) < 0.5: + break + # append results to lists TList = np.append(TList, [ode.t], axis=0) fieldList = np.append(fieldList, [ode.y], axis=0) potentialEffList = np.append(potentialEffList, [potentialEffT], axis=0) + dVdTList = np.append(dVdTList, [[dVdT]], axis=0) + d2VdT2List = np.append(d2VdT2List, [[D2VDT2]], axis=0) + dphidTList = np.append(dphidTList, + [dphidT], + axis=0) if direction == 0: # populating results array TFullList = TList fieldFullList = fieldList potentialEffFullList = potentialEffList + dVdTFullList = dVdTList + d2VdT2FullList = d2VdT2List + dphidTFullList = dphidTList # making new empty array for downwards integration TList = np.empty(0, dtype=float) fieldList = np.empty((0, phase0.numFields()), dtype=float) potentialEffList = np.empty((0, 1), dtype=float) + dVdTList = np.empty((0, 1), dtype=float) + d2VdT2List = np.empty((0, 1), dtype=float) + dphidTList = np.empty((0, phase0.numFields()), dtype=float) else: if len(TList) > 1: # combining up and down integrations @@ -410,6 +463,13 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: potentialEffFullList = np.append( np.flip(potentialEffList, axis=0), potentialEffFullList, axis=0 ) + dVdTFullList = np.append( + np.flip(dVdTList, axis=0), dVdTFullList, axis=0) + d2VdT2FullList = np.append( + np.flip(d2VdT2List, axis=0), d2VdT2FullList, axis=0) + dphidTFullList = np.append( + np.flip(dphidTList, axis=0), dphidTFullList, axis=0 + ) elif len(TFullList) <= 1: # Both up and down lists are too short raise RuntimeError("Failed to trace phase") @@ -439,6 +499,82 @@ def spinodalEvent(temperature: float, field: np.ndarray) -> float: Try decreasing temperatureVariationScale.""" ) + # Compute the second derivative of the field by finite differences + d2phidT2 = np.zeros_like(dphidTFullList) + d2phidT2[1:-1] = ((dphidTFullList[2:] - dphidTFullList[:-2]) / + (TFullList[2:] - TFullList[:-2])[:,None]) + d2phidT2[0] = ((dphidTFullList[1] - dphidTFullList[0]) / + (TFullList[1] - TFullList[0])) + d2phidT2[-1] = ((dphidTFullList[-1] - dphidTFullList[-2]) / + (TFullList[-1] - TFullList[-2])) + # Now to construct the interpolation result = np.concatenate((fieldFullList, potentialEffFullList), axis=1) - self.newInterpolationTableFromValues(TFullList, result) + deriv1 = np.concatenate((dphidTFullList, dVdTFullList), axis=1) + deriv2 = np.concatenate((d2phidT2, d2VdT2FullList), axis=1) + self.newInterpolationTableFromValues(TFullList, result, [deriv1, deriv2], + interpolationDegree) + + def constructInterpolationFromArray( + self, + freeEnergyArrays: "FreeEnergyArrays", + dT: float, + ) -> None: + """ + Constructs the interpolation table directly from arrays of temperatures, + field values, and potential values, bypassing the phase tracing process. + + Parameters + ---------- + freeEnergyArrays : FreeEnergyArrays + Object containing arrays of the temperature, minimum and value of the free energy. + dT : float + Small step in temperature used in derivatives, used here to ensure endpoints of temperature range not exceeded when taking derivatives later. + """ + if freeEnergyArrays.allowedDiscrepancy is None: + freeEnergyArrays.allowedDiscrepancy = self.effectivePotential.effectivePotentialError + + freeEnergyList = freeEnergyArrays.freeEnergyList + + # Check if the loaded value is consistent with the effective potential + discrepancies = abs( + freeEnergyList.veffValue - self.effectivePotential.evaluate( + freeEnergyList.fieldsAtMinimum, freeEnergyArrays.temperatures + ) + ) + # Norm includes neighbouring points to avoid division by zero + discrepancyNorm = np.concatenate( + ([0.5 * (abs(freeEnergyList.veffValue[0]) + abs(freeEnergyList.veffValue[1]))], + 0.5 * (abs(freeEnergyList.veffValue[:-1]) + abs(freeEnergyList.veffValue[1:]))) + ) + maxDiscrepancy = max(discrepancies / discrepancyNorm) + if (maxDiscrepancy > freeEnergyArrays.allowedDiscrepancy): + raise WallGoError( + f"The loaded phase disagrees with the effective potential at {maxDiscrepancy:g}, higher than the required precision of {freeEnergyArrays.allowedDiscrepancy:g}." + ) + + # Check that the provided array has sufficiently small temperature steps + maxTemperatureStep = max(abs( + freeEnergyArrays.temperatures[:-1]-freeEnergyArrays.temperatures[1:] + )) + + if maxTemperatureStep > dT: + logging.warning( + "The maximum temperature step size %g is larger than the maximum step size %g. " + "This may lead to inaccurate interpolation.", + maxTemperatureStep, + dT, + ) + + self.minPossibleTemperature[0] = min(freeEnergyArrays.temperatures) + 2 * dT + self.maxPossibleTemperature[0] = max(freeEnergyArrays.temperatures) - 2 * dT + + # Concatenate field values and potential into a single array (N, nFields + 1) + resultArray = np.concatenate( + (freeEnergyList.fieldsAtMinimum, freeEnergyList.veffValue[:, np.newaxis]), axis=1 + ) + + # Construct the interpolation table + self.newInterpolationTableFromValues( + freeEnergyArrays.temperatures, resultArray + ) diff --git a/src/WallGo/hydrodynamicsTemplateModel.py b/src/WallGo/hydrodynamicsTemplateModel.py index 307faa34..53002262 100644 --- a/src/WallGo/hydrodynamicsTemplateModel.py +++ b/src/WallGo/hydrodynamicsTemplateModel.py @@ -99,6 +99,7 @@ def __init__( self.vMin = self.minVelocity() self.epsilon = self.wN*(1/self.mu-(1-3*self.alN)/self.nu) + def findJouguetVelocity(self, alN: float | None = None) -> float: r""" Finds the Jouguet velocity, corresponding to the phase transition strength diff --git a/src/WallGo/interpolatableFunction.py b/src/WallGo/interpolatableFunction.py index 4643fab0..cc7367c0 100644 --- a/src/WallGo/interpolatableFunction.py +++ b/src/WallGo/interpolatableFunction.py @@ -7,8 +7,7 @@ from typing import Callable, Tuple import logging import numpy as np -from scipy.interpolate import CubicSpline - +from scipy.interpolate import make_interp_spline, BSpline from . import helpers inputType = list[float] | np.ndarray @@ -95,7 +94,7 @@ def __init__( assert returnValueCount >= 1 self._RETURN_VALUE_COUNT = returnValueCount # pylint: disable=invalid-name - self._interpolatedFunction: CubicSpline + self._interpolatedFunction: BSpline ## Will hold list of interpolated derivatives, 1st and 2nd derivatives only self._interpolatedDerivatives: list[Callable] @@ -242,6 +241,8 @@ def newInterpolationTableFromValues( self, x: inputType, fx: outputType, + derivatives: list[outputType] | None = None, + splineDegree: int = 3 ) -> None: """ Like initializeInterpolationTable but takes in precomputed function values 'fx' @@ -252,9 +253,11 @@ def newInterpolationTableFromValues( Points where the function was evaluated. fx : list[float | np.ndarray] or np.ndarray Value of the function at x. - + derivatives : list[outputType] | None + List containing the values of each derivative of the function at x. If None, + computes the derivatives from the interpolated spline. """ - self._interpolate(x, fx) + self._interpolate(x, fx, derivatives, splineDegree) def scheduleForInterpolation(self, x: inputType, fx: outputType) -> None: """ @@ -504,7 +507,9 @@ def derivative( """ x = np.asanyarray(x) - if not bUseInterpolation or not self.hasInterpolation() or order > 2: + if (not bUseInterpolation or + not self.hasInterpolation() or + order > len(self._interpolatedDerivatives)): return helpers.derivative(self._evaluateDirectly, x, n=order) # Use interpolated values whenever possible @@ -555,6 +560,8 @@ def _interpolate( self, x: inputType, fx: outputType, + derivatives: list[outputType] | None = None, + splineDegree: int = 3, ) -> None: """Does the actual interpolation and sets some internal values. Input x needs to be 1D, and input fx needs to be at most 2D. @@ -567,19 +574,21 @@ def _interpolate( ) ## Can't specify different extrapolation methods for x > xmax, x < xmin in - ## CubicSpline! This logic is handled manually in __call__() + ## Spline! This logic is handled manually in __call__() bShouldExtrapolate = EExtrapolationType.FUNCTION in ( self.extrapolationTypeLower, self.extrapolationTypeUpper, ) ## Explicitly drop non-numerics - xFiltered, fxFiltered = self._dropBadPoints(x, fx) + xFiltered, fxFiltered, derivativesFiltered = self._dropBadPoints(x, fx, + derivatives) ## This works even if f(x) is vector valued - self._interpolatedFunction = CubicSpline( - xFiltered, fxFiltered, extrapolate=bShouldExtrapolate, axis=0 + self._interpolatedFunction = make_interp_spline( + xFiltered, fxFiltered, k=splineDegree, axis=0 ) + self._interpolatedFunction.extrapolate = bShouldExtrapolate self._rangeMin = np.min(xFiltered) self._rangeMax = np.max(xFiltered) @@ -589,32 +598,54 @@ def _interpolate( """Store a cubic spline for the 1st and 2nd derivatives into a list. We do not attempt to spline the higher derivatives as they are not guaranteed to be continuous.""" - self._interpolatedDerivatives = [ - self._interpolatedFunction.derivative(1), - self._interpolatedFunction.derivative(2), - ] + if derivatives is None or len(derivatives) == 0: + self._interpolatedDerivatives = [ + self._interpolatedFunction.derivative(1), + self._interpolatedFunction.derivative(2), + ] + else: + self._interpolatedDerivatives = [] + for d in derivativesFiltered: + self._interpolatedDerivatives.append(make_interp_spline( + xFiltered, d, k=splineDegree, axis=0 + )) + self._interpolatedDerivatives[-1].extrapolate = bShouldExtrapolate + if len(self._interpolatedDerivatives) == 1: + self._interpolatedDerivatives.append( + self._interpolatedDerivatives[0].derivative(1)) @staticmethod def _dropBadPoints( x: np.ndarray, fx: np.ndarray, - ) -> tuple[np.ndarray, np.ndarray]: + derivatives: list[outputType] | None = None, + ) -> tuple[np.ndarray, np.ndarray, list[outputType] | None]: """ Removes non-numerical (x, fx) pairs. For 2D fx the check is applied row-wise. Input x needs to be 1D, and input fx needs to be at most 2D. Output is same shape as input. """ + if derivatives is None: + derivativesValid = None + else: + derivativesValid = [] if fx.ndim > 1: validIndices = np.all(np.isfinite(fx), axis=1) fxValid = fx[validIndices] + if derivatives is not None: + for d in derivatives: + derivativesValid.append(d[validIndices]) else: ## fx is 1D array validIndices = np.all(np.isfinite(fx)) fxValid = np.ravel(fx[validIndices]) + if derivatives is not None: + for d in derivatives: + derivativesValid.append(np.ravel(d[validIndices])) xValid = np.ravel(x[validIndices]) - return xValid, fxValid + return xValid, fxValid, derivativesValid def _adaptiveInterpolationUpdate(self) -> None: """ @@ -757,8 +788,11 @@ def readInterpolationTable(self, fileToRead: str) -> None: self._validateInterpolationTable((self._rangeMax - self._rangeMin) / 2.55) logging.debug( - f"{selfName}: Succesfully read interpolation table from file. " - "Range [{self._rangeMin}, {self._rangeMax}]" + "%s: Succesfully read interpolation table from file. " + "Range [%g, %g]", + selfName, + self._rangeMin, + self._rangeMax, ) except IOError as ioError: diff --git a/src/WallGo/manager.py b/src/WallGo/manager.py index dc27d8de..a2d735a4 100644 --- a/src/WallGo/manager.py +++ b/src/WallGo/manager.py @@ -10,7 +10,7 @@ # WallGo imports import WallGo -from .boltzmann import BoltzmannSolver +from .boltzmann import BoltzmannSolver, ETruncationOption from .containers import PhaseInfo from .equationOfMotion import EOM from .exceptions import WallGoError, WallGoPhaseValidationError @@ -68,7 +68,7 @@ class WallSolver: 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. @@ -115,12 +115,14 @@ def setVerbosity(self, verbosityLevel: int) -> None: shown. """ - logging.basicConfig(format='%(message)s', level=verbosityLevel, force=True) + logging.basicConfig(format="%(message)s", level=verbosityLevel, force=True) 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 @@ -158,7 +160,10 @@ def setupThermodynamicsHydrodynamics( # Checks that phase input makes sense with the user-specified Veff self.validatePhaseInput(phaseInfo) - self.initTemperatureRange() + self.initTemperatureRange( + freeEnergyArraysHighT=freeEnergyArraysHighT, + freeEnergyArraysLowT=freeEnergyArraysLowT, + ) ## Should we write these to a result struct? logging.info("Temperature ranges:") @@ -233,13 +238,13 @@ def validatePhaseInput(self, phaseInput: PhaseInfo) -> None: phaseLocation1, effPotValue1, ) = self.model.getEffectivePotential().findLocalMinimum( - phaseInput.phaseLocation1, T + phaseInput.phaseLocation1, T, method='Nelder-Mead' ) ( phaseLocation2, effPotValue2, ) = self.model.getEffectivePotential().findLocalMinimum( - phaseInput.phaseLocation2, T + phaseInput.phaseLocation2, T, method='Nelder-Mead' ) logging.info(f"Found phase 1: phi = {phaseLocation1}, Veff(phi) = {effPotValue1}") @@ -277,11 +282,24 @@ def validatePhaseInput(self, phaseInput: PhaseInfo) -> None: self.phasesAtTn = foundPhaseInfo - def initTemperatureRange(self) -> None: + def initTemperatureRange( + self, + freeEnergyArraysHighT: WallGo.FreeEnergyArrays = None, + freeEnergyArraysLowT: WallGo.FreeEnergyArrays = None, + ) -> None: """ Determine the relevant temperature range and trace the phases over this range. Interpolate the free energy in both phases and store in internal 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 @@ -308,7 +326,9 @@ def initTemperatureRange(self) -> None: # 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()}") + logging.info( + f"vwLTE in the template model: {hydrodynamicsTemplate.findvwLTE()}" + ) except WallGoError as error: # Throw new error with more info @@ -323,6 +343,56 @@ def initTemperatureRange(self) -> None: "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 @@ -341,6 +411,7 @@ def initTemperatureRange(self) -> 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. @@ -348,7 +419,6 @@ def initTemperatureRange(self) -> None: 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. @@ -360,24 +430,23 @@ def initTemperatureRange(self) -> None: TMinLowT = TLowTMinTemplate * self.config.configThermodynamics.tmin TMaxLowT = TLowTMaxTemplate * self.config.configThermodynamics.tmax - # Interpolate phases and check that they remain stable in this range - fHighT = self.thermodynamics.freeEnergyHigh - fLowT = self.thermodynamics.freeEnergyLow - - fHighT.tracePhase( - TMinHighT, - TMaxHighT, - dT, - rTol=phaseTracerTol, - phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep, - ) - fLowT.tracePhase( - TMinLowT, - TMaxLowT, - dT, - rTol=phaseTracerTol, - phaseTracerFirstStep=self.config.configThermodynamics.phaseTracerFirstStep, - ) + # 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, + ) def setPathToCollisionData(self, directoryPath: pathlib.Path) -> None: """ @@ -418,7 +487,7 @@ def solveWall( ) -> 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 @@ -529,13 +598,17 @@ def setupWallSolver(self, wallSolverSettings: WallSolverSettings) -> WallSolver: # 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) @@ -591,13 +664,17 @@ def buildGrid( 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 + tailLength = ( + max( + meanFreePathScale, + 0.5 * wallThicknessIni * (1.0 + 3.0 * smoothing) / ratioPointsWall, + ) + / Tnucl + ) if gridN % 2 == 0: raise ValueError( @@ -617,7 +694,10 @@ def buildGrid( ) def buildEOM( - self, grid: Grid3Scales, boltzmannSolver: BoltzmannSolver, meanFreePathScale: float + self, + grid: Grid3Scales, + boltzmannSolver: BoltzmannSolver, + meanFreePathScale: float, ) -> EOM: r""" Constructs an :py:class:`EOM` object using internal state from the :py:class:`WallGoManager`, @@ -649,7 +729,7 @@ def buildEOM( wallThicknessBounds = self.config.configEOM.wallThicknessBounds wallOffsetBounds = self.config.configEOM.wallOffsetBounds - + Tnucl = self.phasesAtTn.temperature return EOM( @@ -658,7 +738,7 @@ def buildEOM( self.hydrodynamics, grid, numberOfFields, - meanFreePathScale / Tnucl, # We divide by Tnucl to get physical units + meanFreePathScale / Tnucl, # We divide by Tnucl to get physical units wallThicknessBounds, wallOffsetBounds, includeOffEq=True, @@ -668,4 +748,3 @@ def buildEOM( maxIterations=maxIterations, pressRelErrTol=pressRelErrTol, ) - diff --git a/src/WallGo/polynomial.py b/src/WallGo/polynomial.py index fa645c2e..c1c51272 100644 --- a/src/WallGo/polynomial.py +++ b/src/WallGo/polynomial.py @@ -4,6 +4,7 @@ from __future__ import annotations import typing +import logging import numpy as np import numpy.typing as npt from scipy.special import eval_chebyt, eval_chebyu @@ -947,3 +948,90 @@ def _isBroadcastable(self, array1: np.ndarray, array2: np.ndarray) -> bool: else: return False return True + + +class SpectralConvergenceInfo: + """ + Carries information about the convergence of a polynomial expansion. + + Assumes input is a 1d array of coefficients in the Chebyshev basis. Fits a slope to the logarithm of the absolute value of these coefficients. Also, finds the average value of the index, treating the coefficients as a + histogram. + """ + coefficients: np.ndarray + """Coefficients of expansion in the Chebyshev basis, must be 1d.""" + + weightPower: int = 0 + r"""Additional powers of :math:`n` to weight by in assessing convergence. + Default is zero.""" + + offset: int = 0 + r"""Offest in :math:`n`. Default is zero.""" + + apparentConvergence: bool = False + """True if spectral expansion appears to be converging, False otherwise.""" + + spectralPeak: int = 0 + """Positions of the peak of the spectral expansion.""" + + spectralExponent: float = 0.0 + r"""Exponent :math:`\sigma` of :math:`A e^{\sigma n}` fit to spectral expansion.""" + + def __init__( + self, coefficients: np.ndarray, weightPower: int=0, offset: int=0 + ) -> None: + """Initialise given an array.""" + assert len(coefficients.shape) == 1, "SpectralConvergenceInfo requires a 1d array" + self.coefficients = coefficients[:] + self.weightPower = weightPower + self.offset = offset + self._checkSpectralConvergence() + + def _checkSpectralConvergence(self) -> None: + """ + Check for spectral convergence, performing fits and looking at the + position of the maximum. + """ + nCoeff = len(self.coefficients) + weight = (1 + self.offset + np.arange(nCoeff)) ** self.weightPower + absCoefficients = abs(self.coefficients) + if nCoeff < 2: + logging.warning("Spectral convergence tests not valid for n < 2.") + return + if nCoeff < 3: + strict = False + else: + strict = True + + # fit slope to the log of the coefficients + if strict: + # enough points to fit slopes with errors + fit = np.polyfit( + np.arange(nCoeff) + self.offset, + np.log(absCoefficients), + # np.log(absCoefficients * weight), + 1, + cov=True, + ) + self.spectralExponent = fit[0][0] + # self.apparentConvergence = fit[0][0] < - np.sqrt(fit[1][0, 0]) + self.apparentConvergence = fit[0][0] < 0 + else: + # not enough points to get sensible errors + fit = np.polyfit( + np.arange(nCoeff) + self.offset, + np.log(absCoefficients), + # np.log(absCoefficients * weight), + 1, + cov=False, + ) + self.spectralExponent = fit[0] + self.apparentConvergence = fit[0] < 0 + + # Index weighted by absolute value in array + self.spectralPeak = int( + np.sum((np.arange(nCoeff) + self.offset) * weight * absCoefficients) / np.sum(absCoefficients * weight) + ) + + # Alternative convergence condition + # self.apparentConvergence = self.spectralPeak < self.offset + nCoeff // 2 + \ No newline at end of file diff --git a/src/WallGo/results.py b/src/WallGo/results.py index f406ccf5..51adfc97 100644 --- a/src/WallGo/results.py +++ b/src/WallGo/results.py @@ -5,6 +5,7 @@ from dataclasses import dataclass from enum import Enum import numpy as np +from typing import Tuple from .fields import Fields from .containers import BoltzmannBackground, BoltzmannDeltas, WallParams @@ -27,24 +28,31 @@ class BoltzmannResults: r"""Estimated relative error in :math:`\delta f` due to truncation of spectral expansion.""" + truncatedTail: tuple[bool, bool, bool] + """True if tail 1/3 of spectral expansion was truncated in each direction, False otherwise.""" + + spectralPeaks: tuple[int, int, int] + r"""Indices of the peaks in the spectral expansion in each direction.""" + # These two criteria are to evaluate the validity of the linearization of the # Boltzmann equation. The arrays contain one element for each out-of-equilibrium # particle. To be valid, at least one criterion must be small for each particle. - linearizationCriterion1: np.ndarray + linearizationCriterion1: float = 0 r"""Ratio of out-of-equilibrium and equilibrium pressures, - :math:`|P[\delta f]| / |P[f_\text{eq}]|`. One element for each - out-of-equilibrium particle.""" + :math:`|P[\delta f]| / |P[f_\text{eq}]|`. Default is 0.""" - linearizationCriterion2: np.ndarray - r"""Ratio of collision and Liouville operators in Boltzmann equation, - :math:`|\mathcal{C}[\delta f]|/ |\mathcal{L}[\delta f]|`. One element for each - out-of-equilibrium particle.""" + linearizationCriterion2: float = 0 + r"""Ratio of the first-order correction due to nonlinearities and total pressure + computed by WallGo, :math:`|P[\delta f_2]| / |P[f_\text{eq}+\delta f]|`. + Default is 0.""" def __mul__(self, number: float) -> "BoltzmannResults": return BoltzmannResults( deltaF=number * self.deltaF, Deltas=number * self.Deltas, truncationError=abs(number) * self.truncationError, + truncatedTail=self.truncatedTail, + spectralPeaks=self.spectralPeaks, linearizationCriterion1=abs(number) * self.linearizationCriterion1, linearizationCriterion2=self.linearizationCriterion2, ) @@ -54,6 +62,8 @@ def __rmul__(self, number: float) -> "BoltzmannResults": deltaF=number * self.deltaF, Deltas=number * self.Deltas, truncationError=abs(number) * self.truncationError, + truncatedTail=self.truncatedTail, + spectralPeaks=self.spectralPeaks, linearizationCriterion1=abs(number) * self.linearizationCriterion1, linearizationCriterion2=self.linearizationCriterion2, ) @@ -63,6 +73,8 @@ def __add__(self, other: "BoltzmannResults") -> "BoltzmannResults": deltaF=other.deltaF + self.deltaF, Deltas=other.Deltas + self.Deltas, truncationError=other.truncationError + self.truncationError, + truncatedTail=self.truncatedTail, + spectralPeaks=self.spectralPeaks, linearizationCriterion1=other.linearizationCriterion1 + self.linearizationCriterion1, linearizationCriterion2=other.linearizationCriterion2 @@ -172,15 +184,29 @@ class WallGoResults: temperatureProfile: np.ndarray r"""Temperarture profile as a function of position, :math:`T(\xi)`.""" - linearizationCriterion1: np.ndarray + linearizationCriterion1: float r"""Ratio of out-of-equilibrium and equilibrium pressures, - :math:`|P[\delta f]| / |P[f_\text{eq}]|`. One element for each - out-of-equilibrium particle.""" + :math:`|P[\delta f]| / |P[f_\text{eq}]|`.""" + + eomResidual: np.ndarray + r""" + Residual of the EOM due to the tanh ansatz. There is one element for each scalar + field. It is estimated by the integral + + .. math:: \sqrt{\Delta[\mathrm{EOM}^2]/|\mathrm{EOM}^2|} + + with + + .. math:: \Delta[\mathrm{EOM}^2]=\int\! dz\, (-\partial_z^2 \phi+ \partial V_{\mathrm{eq}}/ \partial \phi+ \partial V_{\mathrm{out}}/ \partial \phi )^2 + + and + + .. math:: |\mathrm{EOM}^2|=\int\! dz\, [(\partial_z^2 \phi)^2+ (\partial V_{\mathrm{eq}}/ \partial \phi)^2+ (\partial V_{\mathrm{out}}/ \partial \phi)^2]. + """ - linearizationCriterion2: np.ndarray - r"""Ratio of collision and Liouville operators in Boltzmann equation, - :math:`|\mathcal{C}[\delta f]|/ |\mathcal{L}[\delta f]|`. One element for each - out-of-equilibrium particle.""" + linearizationCriterion2: float + r"""Ratio of the first-order correction due to nonlinearities and total pressure + computed by WallGo, :math:`|P[\delta f_2]| / |P[f_\text{eq}+\delta f]|`.""" deltaF: np.ndarray r"""Deviation of probability density function from equilibrium, @@ -204,6 +230,12 @@ class WallGoResults: :math:`\mathcal{E}_\text{pl}^{n_\mathcal{E}}\mathcal{P}_\text{pl}^{n_\mathcal{P}}\delta f`, using finite differences instead of spectral expansion.""" + violationOfEMConservation: Tuple[float,float] + """ + RMS along the grid of the violation of the conservation of the components T30 and T33 of + the energy-momentum tensor. + """ + solutionType: ESolutionType """ Describes the type of solution obtained. Must be a ESolutionType object. The @@ -281,6 +313,17 @@ def setFiniteDifferenceBoltzmannResults( self.deltaFFiniteDifference = boltzmannResults.deltaF self.DeltasFiniteDifference = boltzmannResults.Deltas + def setViolationOfEMConservation( + self, violationOfEMConservation: Tuple[float, float] + ) -> None: + """ + Set the violation of energy-momentum conservation results + """ + assert ( + len(violationOfEMConservation) == 2 + ), "WallGoResults Error: violationOfEMConservation must be a tuple of two floats." + self.violationOfEMConservation = violationOfEMConservation + def setSuccessState( self, success: bool, @@ -297,4 +340,4 @@ def setSuccessState( ) self.success = success self.message = message - self.solutionType = solutionType + self.solutionType = solutionType \ No newline at end of file diff --git a/tests/BenchmarkPoint.py b/tests/BenchmarkPoint.py index 687c5a6f..ef036b16 100644 --- a/tests/BenchmarkPoint.py +++ b/tests/BenchmarkPoint.py @@ -1,6 +1,8 @@ import WallGo.genericModel +from typing import Any + ## Collect input params + other benchmark-specific data for various things in one place. class BenchmarkPoint: @@ -10,7 +12,10 @@ class BenchmarkPoint: phaseInfo: dict[str, float] ## This is WallGo internal config info that we may want to fix on a per-benchmark basis. IE. temperature interpolation ranges - config: dict[str, float] + # BREAKING CHANGE: The type annotation for `config` was changed from `dict[str, float]` to `dict[str, Any]` + # to allow for more flexible configuration values (e.g., nested dictionaries for derivative settings). + # This broadens the accepted types and may break code that expects all config values to be floats. + config: dict[str, Any] ## Expected results for the benchmark point expectedResults: dict[str, float] @@ -18,14 +23,14 @@ class BenchmarkPoint: def __init__( self, inputParams: dict[str, float], - phaseInfo: dict[str, float] = {}, - config: dict[str, float] = {}, - expectedResults: dict[str, float] = {}, + phaseInfo: dict[str, float] | None = None, + config: dict[str, Any] | None = None, + expectedResults: dict[str, float] | None = None, ): self.inputParams = inputParams - self.phaseInfo = phaseInfo - self.config = config - self.expectedResults = expectedResults + self.phaseInfo = phaseInfo or {} + self.config = config or {} + self.expectedResults = expectedResults or {} class BenchmarkModel: @@ -36,6 +41,26 @@ class BenchmarkModel: def __init__(self, model: WallGo.GenericModel, benchmarkPoint: BenchmarkPoint): self.model = model - self.model.getEffectivePotential().configureDerivatives(WallGo.VeffDerivativeSettings(1.0, 1.0)) + + # Apply derivative settings if specified in benchmark configuration + derivative_settings = benchmarkPoint.config.get("derivativeSettings") + if derivative_settings is not None: + if isinstance(derivative_settings, dict): + temp_scale = derivative_settings.get("temperatureVariationScale", 1.0) + field_scale = derivative_settings.get("fieldValueVariationScale", 1.0) + self.model.getEffectivePotential().configureDerivatives( + WallGo.VeffDerivativeSettings(temp_scale, field_scale) + ) + else: + # Use default settings if configuration is malformed + self.model.getEffectivePotential().configureDerivatives( + WallGo.VeffDerivativeSettings(1.0, 1.0) + ) + else: + # Use default settings if not specified + self.model.getEffectivePotential().configureDerivatives( + WallGo.VeffDerivativeSettings(1.0, 1.0) + ) + self.model.getEffectivePotential().effectivePotentialError = 1e-15 self.benchmarkPoint = benchmarkPoint diff --git a/tests/Benchmarks/SingletSM_Z2/conftest.py b/tests/Benchmarks/SingletSM_Z2/conftest.py index 7f2bb11b..42efa531 100644 --- a/tests/Benchmarks/SingletSM_Z2/conftest.py +++ b/tests/Benchmarks/SingletSM_Z2/conftest.py @@ -123,8 +123,11 @@ def singletBenchmarkThermo_interpolate( """ Then manually interpolate """ TMin, TMax, dT = BM.config["interpolateTemperatureRange"] - thermo.freeEnergyHigh.tracePhase(TMin, TMax, dT) - thermo.freeEnergyLow.tracePhase(TMin, TMax, dT) + # To meet the high accuracy requirement of this test, we set the interpolation order + # to 3. We do not recommend to do this in general, as it can lead to unphysical + #features in the speed of sound. + thermo.freeEnergyHigh.tracePhase(TMin, TMax, dT, interpolationDegree=3) + thermo.freeEnergyLow.tracePhase(TMin, TMax, dT, interpolationDegree=3) thermo.setExtrapolate() @@ -170,8 +173,11 @@ def singletSimpleBenchmarkFreeEnergy( dT = 0.1 BM.config["interpolateTemperatureRange"] = TMin, TMax, dT - freeEnergy1.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False) - freeEnergy2.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False) + # To meet the high accuracy requirement of this test, we set the interpolation order + # to 3. We do not recommend to do this in general, as it can lead to unphysical + #features in the speed of sound. + freeEnergy1.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False, interpolationDegree=3) + freeEnergy2.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False, interpolationDegree=3) yield freeEnergy1, freeEnergy2, BM @@ -207,8 +213,11 @@ def singletSimpleBenchmarkThermodynamics( thermo.freeEnergyHigh.disableAdaptiveInterpolation() thermo.freeEnergyLow.disableAdaptiveInterpolation() - thermo.freeEnergyHigh.tracePhase(TMin, TMax, dT) - thermo.freeEnergyLow.tracePhase(TMin, TMax, dT) + # To meet the high accuracy requirement of this test, we set the interpolation order + # to 3. We do not recommend to do this in general, as it can lead to unphysical + #features in the speed of sound. + thermo.freeEnergyHigh.tracePhase(TMin, TMax, dT, interpolationDegree=3) + thermo.freeEnergyLow.tracePhase(TMin, TMax, dT, interpolationDegree=3) thermo.setExtrapolate() diff --git a/tests/Benchmarks/SingletSM_Z2/test_FreeEnergy.py b/tests/Benchmarks/SingletSM_Z2/test_FreeEnergy.py index e96e0a65..332a31b9 100644 --- a/tests/Benchmarks/SingletSM_Z2/test_FreeEnergy.py +++ b/tests/Benchmarks/SingletSM_Z2/test_FreeEnergy.py @@ -2,7 +2,7 @@ import numpy as np from typing import Tuple -from tests.BenchmarkPoint import BenchmarkPoint +from tests.BenchmarkPoint import BenchmarkPoint, BenchmarkModel import WallGo @@ -11,7 +11,7 @@ def test_freeEnergy_singletSimple( singletSimpleBenchmarkFreeEnergy: Tuple[WallGo.FreeEnergy, WallGo.FreeEnergy, BenchmarkPoint], T: float, -): +) -> None: """ Testing numerics of FreeEnergy """ @@ -47,3 +47,77 @@ def test_freeEnergy_singletSimple( assert vExact == pytest.approx(v, rel=rTol) assert 0 == pytest.approx(x, abs=aTol) assert f0 + VvExact == pytest.approx(veffValue, rel=rTol) + + +def test_freeEnergy_singletSimple_passingArrays( + singletSimpleBenchmarkModel: BenchmarkModel, + singletSimpleBenchmarkFreeEnergy: Tuple[WallGo.FreeEnergy, WallGo.FreeEnergy, BenchmarkPoint], +) -> None: + """ + Testing building FreeEnergy from passing arrays + """ + freeEnergy1, freeEnergy2, BM = singletSimpleBenchmarkFreeEnergy + + # temperature range + temperatureRange = np.linspace(90, 110, num=50) + nT = len(temperatureRange) + + vList = np.zeros((nT, 2)) + vVeffList = np.zeros(nT) + xList = np.zeros((nT, 2)) + xVeffList = np.zeros(nT) + + # tolerance + tol = 1e-15 + + for iT, T in enumerate(temperatureRange): + # exact results + thermalParameters = freeEnergy1.effectivePotential.getThermalParameters(T) + f0 = -107.75 * np.pi ** 2 / 90 * T ** 4 + + vExact = np.sqrt(-thermalParameters["muHsq"] / thermalParameters["lHH"]) + VvExact = -0.25 * thermalParameters["muHsq"] ** 2 / thermalParameters["lHH"] + vList[iT, 0] = vExact + vVeffList[iT] = f0 + VvExact + + xExact = np.sqrt(-thermalParameters["muSsq"] / thermalParameters["lSS"]) + VxExact = -0.25 * thermalParameters["muSsq"] ** 2 / thermalParameters["lSS"] + + xList[iT, 1] = xExact + xVeffList[iT] = f0 + VxExact + + freeEnergyHighT = WallGo.FreeEnergyArrays( + temperatures=temperatureRange, + minimumList=xList, + potentialEffList=xVeffList, + allowedDiscrepancy=tol, + ) + + freeEnergyLowT = WallGo.FreeEnergyArrays( + temperatures=temperatureRange, + minimumList=vList, + potentialEffList=vVeffList, + allowedDiscrepancy=tol, + ) + + # free energies for both phases + freeEnergy1 = WallGo.FreeEnergy( + singletSimpleBenchmarkModel.model.getEffectivePotential(), + temperatureRange[10], + WallGo.Fields(vList[10]), + ) + freeEnergy2 = WallGo.FreeEnergy( + singletSimpleBenchmarkModel.model.getEffectivePotential(), + temperatureRange[10], + WallGo.Fields(xList[10]), + ) + + freeEnergy1.constructInterpolationFromArray( + freeEnergyArrays=freeEnergyLowT, + dT=abs(temperatureRange[1]-temperatureRange[0]) + ) + + freeEnergy2.constructInterpolationFromArray( + freeEnergyArrays=freeEnergyHighT, + dT=abs(temperatureRange[1]-temperatureRange[0]) + ) diff --git a/tests/Benchmarks/StandardModel/Benchmarks_SM.py b/tests/Benchmarks/StandardModel/Benchmarks_SM.py new file mode 100644 index 00000000..3ea49d00 --- /dev/null +++ b/tests/Benchmarks/StandardModel/Benchmarks_SM.py @@ -0,0 +1,222 @@ +import WallGo.fields + +from tests.BenchmarkPoint import BenchmarkPoint + +BM1 = BenchmarkPoint( + inputParams={ + "v0": 246.0, + "mW": 80.4, + "mZ": 91.2, + "mt": 174.0, + "g3": 1.2279920495357861, + "mH": 0.0, + }, + phaseInfo={ + "Tn":57.1958, + ## Guesses for phase locations + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([55.]), + }, + config={ + ## Give TMin, TMax, dT as a tuple + "interpolateTemperatureRangeHighTPhase": ( + 46, + 58., + 0.01, + ), ##Set by hand such that the lowest temperature is (slightly) below the actual minimum temperature + "interpolateTemperatureRangeLowTPhase": ( + 56., + 58., + 0.01, + ), ##Set by hand such that the highest temperature is (slightly) above the actual maximum temperature + "derivativeSettings": { + "temperatureVariationScale": 0.75, + "fieldValueVariationScale": 50.0, + }, ## Optimized settings for Standard Model + }, + ## Will probs need to adjust these once we decide on what our final implementation of the benchmark model is + expectedResults={ + "Tc": 57.4104, + ## Phase locations at nucleation temperature + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([55.806]), + "minimumTemperaturePhase1": 51.8199, + "maximumTemperaturePhase2": 57.583, + }, +) + +BM2 = BenchmarkPoint( + inputParams={ + "v0": 246.0, + "mW": 80.4, + "mZ": 91.2, + "mt": 174.0, + "g3": 1.2279920495357861, + "mH": 34.0, + }, + phaseInfo={ + "Tn": 70.5793, + ## Guesses for phase locations + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([64.6294]), + }, + config={ + ## Give TMin, TMax, dT as a tuple + "interpolateTemperatureRangeHighTPhase": ( + 60., + 75., + 0.01, + ), ##Set by hand such that the lowest temperature is (slightly) below the actual minimum temperature + "interpolateTemperatureRangeLowTPhase": ( + 66., + 72.5, + 0.01, + ), ##Set by hand such that the highest temperature is (slightly) above the actual maximum temperature + "derivativeSettings": { + "temperatureVariationScale": 0.75, + "fieldValueVariationScale": 50.0, + }, ## Optimized settings for Standard Model + }, + ## Will probs need to adjust these once we decide on what our final implementation of the benchmark model is + expectedResults={ + "Tc": 70.8238, + ## Phase locations at nucleation temperature + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([64.9522]), + "minimumTemperaturePhase1": 63.8464, + "maximumTemperaturePhase2": 71.025, + }, +) + +BM3 = BenchmarkPoint( + inputParams={ + "v0": 246.0, + "mW": 80.4, + "mZ": 91.2, + "mt": 174.0, + "g3": 1.2279920495357861, + "mH": 50.0, + }, + phaseInfo={ + "Tn": 83.4251, + ## Guesses for phase locations + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([67.2538]), + }, + config={ + ## Give TMin, TMax, dT as a tuple + "interpolateTemperatureRangeHighTPhase": ( + 70., + 90., + 0.01, + ), ##Set by hand such that the lowest temperature is (slightly) below the actual minimum temperature + "interpolateTemperatureRangeLowTPhase": ( + 75., + 85., + 0.01, + ), ##Set by hand such that the highest temperature is (slightly) above the actual maximum temperature + "derivativeSettings": { + "temperatureVariationScale": 0.75, + "fieldValueVariationScale": 50.0, + }, ## Optimized settings for Standard Model + }, + ## Will probs need to adjust these once we decide on what our final implementation of the benchmark model is + expectedResults={ + "Tc": 83.668, + ## Phase locations at nucleation temperature + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([67.2538]), + "minimumTemperaturePhase1": 75.291, + "maximumTemperaturePhase2": 83.879, + }, +) + +BM4 = BenchmarkPoint( + inputParams={ + "v0": 246.0, + "mW": 80.4, + "mZ": 91.2, + "mt": 174.0, + "g3": 1.2279920495357861, + "mH": 70.0, + }, + phaseInfo={ + "Tn": 102.344, + ## Guesses for phase locations + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([65.8969]), + }, + config={ + ## Give TMin, TMax, dT as a tuple + "interpolateTemperatureRangeHighTPhase": ( + 95., + 110., + 0.01, + ), ##Set by hand such that the lowest temperature is (slightly) below the actual minimum temperature + "interpolateTemperatureRangeLowTPhase": ( + 90., + 105., + 0.01, + ), ##Set by hand such that the highest temperature is (slightly) above the actual maximum temperature + "derivativeSettings": { + "temperatureVariationScale": 0.75, + "fieldValueVariationScale": 50.0, + }, ## Optimized settings for Standard Model + }, + ## Will probs need to adjust these once we decide on what our final implementation of the benchmark model is + expectedResults={ + "Tc": 102.57, + ## Phase locations at nucleation temperature + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([65.8969]), + "minimumTemperaturePhase1": 91.9013, + "maximumTemperaturePhase2": 102.844, + }, +) + +BM5 = BenchmarkPoint( + inputParams={ + "v0": 246.0, + "mW": 80.4, + "mZ": 91.2, + "mt": 174.0, + "g3": 1.2279920495357861, + "mH": 81.0, + }, + phaseInfo={ + "Tn": 113.575, + ## Guesses for phase locations + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([64.8777]), + }, + config={ + ## Give TMin, TMax, dT as a tuple + "interpolateTemperatureRangeHighTPhase": ( + 105., + 118., + 0.01, + ), ##Set by hand such that the lowest temperature is (slightly) below the actual minimum temperature + "interpolateTemperatureRangeLowTPhase": ( + 105., + 116., + 0.01, + ), ##Set by hand such that the highest temperature is (slightly) above the actual maximum temperature + "derivativeSettings": { + "temperatureVariationScale": 0.75, + "fieldValueVariationScale": 50.0, + }, ## Optimized settings for Standard Model + }, + ## Will probs need to adjust these once we decide on what our final implementation of the benchmark model is + expectedResults={ + "Tc": 113.795, + ## Phase locations at nucleation temperature + "phaseLocation1": WallGo.Fields([0.0]), + "phaseLocation2": WallGo.Fields([64.8777]), + "minimumTemperaturePhase1": 101.602, + "maximumTemperaturePhase2": 114.105, + }, +) + + +## +standardModelBenchmarks = [BM1, BM2, BM3, BM4, BM5] diff --git a/tests/Benchmarks/StandardModel/__init__.py b/tests/Benchmarks/StandardModel/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/Benchmarks/StandardModel/conftest.py b/tests/Benchmarks/StandardModel/conftest.py new file mode 100644 index 00000000..3928613f --- /dev/null +++ b/tests/Benchmarks/StandardModel/conftest.py @@ -0,0 +1,148 @@ +## StandardModel/conftest.py -- Configure singlet model specific tests. +import pytest +from typing import Tuple + +import WallGo + +from tests.BenchmarkPoint import BenchmarkPoint, BenchmarkModel + +from .Benchmarks_SM import standardModelBenchmarks + +from Models.StandardModel.standardModel import ( + StandardModel, +) + + +""" +FreeEnergy interpolations are initialized in the standardModelBenchmarkThermo_interpolate() fixture below. +Interpolations have a huge impact on performance but also affect the results somewhat. +Bottom line is that these tests ARE sensitive to details of the interpolations! +""" + + +"""----- Define fixtures for the standard model. + +NOTE: I'm giving these session scope so that their state is preserved between tests (cleared when pytest finishes). +This is helpful as things like FreeEnergy interpolations are slow, however it does make our tests a bit less transparent. +""" + + +@pytest.fixture(scope="session", params=standardModelBenchmarks, ids=["BM1", "BM2", "BM3", "BM4", "BM5"]) +def standardModelBenchmarkPoint(request: pytest.FixtureRequest) -> BenchmarkPoint: + yield request.param + + +@pytest.fixture(scope="session") +def standardModelBenchmarkModel(standardModelBenchmarkPoint: BenchmarkPoint) -> BenchmarkModel: + inputs = standardModelBenchmarkPoint.inputParams + model = StandardModel() + model.updateModel(inputs) + + yield BenchmarkModel(model, standardModelBenchmarkPoint) + + + +"""----- Fixtures for more complicated things that depend on the model/Veff. +I'm making these return also the original benchmark point so that it's easier to validate results, +eg. read from BenchmarkPoint.expectedResults""" + + +## This constructs thermodynamics without interpolating anything +@pytest.fixture(scope="session") +def standardModelBenchmarkThermo( + standardModelBenchmarkModel: BenchmarkModel, +) -> Tuple[WallGo.Thermodynamics, BenchmarkPoint]: + + BM = standardModelBenchmarkModel.benchmarkPoint + + Tn = BM.phaseInfo["Tn"] + phase1 = BM.expectedResults["phaseLocation1"] + phase2 = BM.expectedResults["phaseLocation2"] + + # I assume phase1 = high-T, phase2 = low-T. Would prefer to drop these labels though, + # so WallGo could safely assume that the transition is always phase1 -> phase2 + thermo = WallGo.Thermodynamics( + standardModelBenchmarkModel.model.getEffectivePotential(), + Tn, + phase2, + phase1, + ) + + thermo.freeEnergyHigh.disableAdaptiveInterpolation() + thermo.freeEnergyLow.disableAdaptiveInterpolation() + + thermo.freeEnergyHigh.minPossibleTemperature = 40.0 + thermo.freeEnergyHigh.maxPossibleTemperature = 60.0 + thermo.freeEnergyLow.minPossibleTemperature = 40.0 + thermo.freeEnergyLow.maxPossibleTemperature = 60.0 + + thermo.setExtrapolate() + + yield thermo, BM + + +## This is like the standardModelBenchmarkThermo fixture but interpolates the FreeEnergy objects over the temperature range specified in our BM input +@pytest.fixture(scope="session") +def standardModelBenchmarkThermo_interpolate( + standardModelBenchmarkModel: BenchmarkModel, +) -> Tuple[WallGo.Thermodynamics, BenchmarkPoint]: + + BM = standardModelBenchmarkModel.benchmarkPoint + + Tn = BM.phaseInfo["Tn"] + phase1 = BM.expectedResults["phaseLocation1"] + phase2 = BM.expectedResults["phaseLocation2"] + + ## I assume phase1 = high-T, phase2 = low-T. Would prefer to drop these labels though, + ## so WallGo could safely assume that the transition is always phase1 -> phase2 + thermo = WallGo.Thermodynamics( + standardModelBenchmarkModel.model.getEffectivePotential(), Tn, phase2, phase1 + ) + + ## Let's turn these off so that things are more transparent + thermo.freeEnergyHigh.disableAdaptiveInterpolation() + thermo.freeEnergyLow.disableAdaptiveInterpolation() + + """ Then manually interpolate """ + TMin, TMax, dT = BM.config["interpolateTemperatureRangeHighTPhase"] + + thermo.freeEnergyHigh.tracePhase(TMin, TMax, dT) + + TMin, TMax, dT = BM.config["interpolateTemperatureRangeLowTPhase"] + thermo.freeEnergyLow.tracePhase(TMin, TMax, dT) + + thermo.setExtrapolate() + + yield thermo, BM + + + +## Test for following minimum +@pytest.fixture(scope="session") +def standardModelBenchmarkFreeEnergy( + standardModelBenchmarkModel: BenchmarkModel, +) -> Tuple[WallGo.FreeEnergy, WallGo.FreeEnergy, BenchmarkPoint]: + + BM = standardModelBenchmarkModel.benchmarkPoint + + Tn = BM.phaseInfo["Tn"] + phase1 = BM.expectedResults["phaseLocation1"] + phase2 = BM.expectedResults["phaseLocation2"] + + # free energies for both phases + freeEnergy1 = WallGo.FreeEnergy( + standardModelBenchmarkModel.model.getEffectivePotential(), Tn, phase1 + ) + freeEnergy2 = WallGo.FreeEnergy( + standardModelBenchmarkModel.model.getEffectivePotential(), Tn, phase2 + ) + + + """ Then manually interpolate """ + TMin, TMax, dT = BM.config["interpolateTemperatureRangeHighTPhase"] + freeEnergy1.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False) + + TMin, TMax, dT = BM.config["interpolateTemperatureRangeLowTPhase"] + freeEnergy2.tracePhase(TMin, TMax, dT, rTol=1e-6, paranoid=False) + + yield freeEnergy1, freeEnergy2, BM \ No newline at end of file diff --git a/tests/Benchmarks/StandardModel/test_Tc.py b/tests/Benchmarks/StandardModel/test_Tc.py new file mode 100644 index 00000000..d6e7b289 --- /dev/null +++ b/tests/Benchmarks/StandardModel/test_Tc.py @@ -0,0 +1,24 @@ +import pytest +import numpy as np +from typing import Tuple + +from tests.BenchmarkPoint import BenchmarkPoint + +import WallGo + +@pytest.mark.slow +def test_standardModelThermodynamicsFindCriticalTemperature( + standardModelBenchmarkThermo_interpolate: Tuple[WallGo.Thermodynamics, BenchmarkPoint], +): + + thermodynamics, BM = standardModelBenchmarkThermo_interpolate + + Tc = thermodynamics.findCriticalTemperature( + dT=0.01, + rTol=1e-8, + paranoid=True, + ) + + expectedTc = BM.expectedResults["Tc"] + + assert Tc == pytest.approx(expectedTc, rel=1e-4) diff --git a/tests/Benchmarks/StandardModel/test_phaseTrace.py b/tests/Benchmarks/StandardModel/test_phaseTrace.py new file mode 100644 index 00000000..434e90db --- /dev/null +++ b/tests/Benchmarks/StandardModel/test_phaseTrace.py @@ -0,0 +1,23 @@ +import pytest +import numpy as np +from typing import Tuple + +from tests.BenchmarkPoint import BenchmarkPoint + +import WallGo + +@pytest.mark.slow +def test_standardModelPhaseTrace( + standardModelBenchmarkFreeEnergy: Tuple[WallGo.FreeEnergy, WallGo.FreeEnergy, BenchmarkPoint], +): + + + _, freeEnergyLowT, BM = standardModelBenchmarkFreeEnergy + + _,_,dT = BM.config["interpolateTemperatureRangeHighTPhase"] + + maxT = freeEnergyLowT.maxPossibleTemperature[0] + 2*dT + expectedMaxT = BM.expectedResults["maximumTemperaturePhase2"] + + + assert maxT == pytest.approx(expectedMaxT, rel=1e-3) \ No newline at end of file diff --git a/tests/test_Boltzmann.py b/tests/test_Boltzmann.py index a4059196..c7fc82cd 100644 --- a/tests/test_Boltzmann.py +++ b/tests/test_Boltzmann.py @@ -38,6 +38,7 @@ def test_Delta00( grid = WallGo.grid.Grid(spatialGridSize, momentumGridSize, 1, 100) collisionPath = dir_path / "TestData/N19" boltzmann = WallGo.BoltzmannSolver(grid, "Cardinal", "Cardinal", "Spectral") + boltzmann.truncationOption = WallGo.ETruncationOption.NONE boltzmann.updateParticleList([particle]) boltzmann.setBackground(bg) @@ -122,3 +123,52 @@ def test_solution( # asserting solution works assert ratio == pytest.approx(0, abs=1e-14) + + +@pytest.mark.parametrize("spatialGridSize, momentumGridSize, slope", [(7, 9, -0.1), (9, 9, -0.1), (11, 9, 0.1)]) +def test_checkSpectralConvergence( + boltzmannTestBackground: WallGo.BoltzmannBackground, + particle: WallGo.Particle, + spatialGridSize: int, + momentumGridSize: int, + slope: float, +) -> None: + """ + Tests that the Boltzmann equation is satisfied by the solution + """ + # setting up objects + grid = WallGo.grid.Grid(spatialGridSize, momentumGridSize, 1, 1) + boltzmann = WallGo.BoltzmannSolver( + grid, + basisM="Chebyshev", + basisN="Chebyshev", + truncationOption=WallGo.ETruncationOption.AUTO, + ) + + # solving Boltzmann equations + deltaFShape = ( + 1, + spatialGridSize - 1, + momentumGridSize - 1, + momentumGridSize - 1, + ) + deltaF = np.zeros(deltaFShape) + for z in range(spatialGridSize - 1): + for pz in range(momentumGridSize - 1): + for pp in range(momentumGridSize - 1): + deltaF[0, z, pz, pp] = np.exp(slope*z - 1e-2*pz - 1e-3*pp) / (1 + pz) / (1 + pp)**2 + + + # checking spectral convergence + deltaFTruncated, truncatedShape, _ = boltzmann.checkSpectralConvergence(deltaF) + + nTruncated = (spatialGridSize - 1) // 3 + expectTruncated = deltaFTruncated[:, -nTruncated:, :, :] + + # asserting truncation + if slope > 0: + assert truncatedShape == (1, spatialGridSize - 1 - nTruncated, momentumGridSize - 1, momentumGridSize - 1) + assert np.allclose(expectTruncated, np.zeros_like(expectTruncated), atol=1e-2) + else: + assert truncatedShape == (1, spatialGridSize - 1, momentumGridSize - 1, momentumGridSize - 1) + assert not np.allclose(expectTruncated, np.zeros_like(expectTruncated)) diff --git a/tests/test_EOM.py b/tests/test_EOM.py new file mode 100644 index 00000000..a5c5d866 --- /dev/null +++ b/tests/test_EOM.py @@ -0,0 +1,324 @@ +import pytest +from dataclasses import dataclass +import numpy as np +from scipy.integrate import odeint +import WallGo + +class QuarticZ2(WallGo.GenericModel): + r""" + Z2 symmetric quartic potential. + + The potential is given by: + V = 1/2 muSq |phi|^2 + 1/4 lam |phi|^4 + + This class inherits from the GenericModel class and implements the necessary + methods for the WallGo package. + """ + + def __init__(self, modelParameters: dict[str, float]): + """ + Initialize the QuarticZ2 model. + """ + + self.modelParameters = modelParameters + + # Initialize internal effective potential + self.effectivePotential = EffectivePotentialQuarticZ2(self) + + # Create a list of particles relevant for the Boltzmann equations + self.defineParticles() + + # ~ GenericModel interface + @property + def fieldCount(self) -> int: + """How many classical background fields""" + return 1 + + def getEffectivePotential(self) -> "EffectivePotentialQuarticZ2": + return self.effectivePotential + + # ~ + + def defineParticles(self) -> None: + """ + Define the particles for the model. + Note that the particle list only needs to contain the + particles that are relevant for the Boltzmann equations. + The particles relevant to the effective potential are + included independently. + + Parameters + ---------- + None + + Returns + ---------- + None + """ + self.clearParticles() + + # === Top quark === + # The msqVacuum function of an out-of-equilibrium particle must take + # a Fields object and return an array of length equal to the number of + # points in fields. + def topMsqVacuum(fields: WallGo.Fields) -> WallGo.Fields: + return 0.5 * fields.getField(0) ** 2 + + # The msqDerivative function of an out-of-equilibrium particle must take + # a Fields object and return an array with the same shape as fields. + def topMsqDerivative(fields: WallGo.Fields) -> WallGo.Fields: + return np.transpose( + [fields.getField(0)] + ) + + topQuark = WallGo.Particle( + "top", + index=0, + msqVacuum=topMsqVacuum, + msqDerivative=topMsqDerivative, + statistics="Fermion", + totalDOFs=12, + ) + self.addParticle(topQuark) + + +# end model + + +class EffectivePotentialQuarticZ2(WallGo.EffectivePotential): + """ + Effective potential for the QuarticZ2 model. + + This class inherits from the EffectivePotential class and provides the + necessary methods for calculating the effective potential. + """ + + # ~ EffectivePotential interface + fieldCount = 1 + """How many classical background fields""" + + effectivePotentialError = 1e-15 + """ + Relative accuracy at which the potential can be computed. + """ + + def __init__(self, owningModel: QuarticZ2) -> None: + """ + Initialize the EffectivePotentialQuarticZ2. + """ + + assert owningModel is not None, "Invalid model passed to Veff" + + self.owner = owningModel + self.modelParameters = self.owner.modelParameters + + # Count particle degrees-of-freedom to facilitate inclusion of + # light particle contributions to ideal gas pressure + self.numBosonDof = 29 + self.numFermionDof = 90 + + def evaluate( + self, fields: WallGo.Fields, temperature: float + ) -> float | np.ndarray: + """ + Evaluate the effective potential. + + Parameters + ---------- + fields: Fields + The field configuration + temperature: float + The temperature + + Returns + ---------- + potentialTotal: complex | np.ndarray + The value of the effective potential + """ + + # phi ~ 1/sqrt(2) (0, v), S ~ x + fields = WallGo.Fields(fields) + v = fields.getField(0) + + muSq = self.modelParameters["muSq"] + lam = self.modelParameters["lam"] + # We include a small cubic term to make the second minimum slightly deeper. + # Otherwise WallGo would not converge. + smallCubic = -1e-4 + + # tree level potential + potentialTree = ( + 0.5 * muSq * v**2 + + 0.25 * lam * v**4 + + smallCubic * v**3 / 3 + ) + + # Include only the T^4 terms in the potential + potentialTotal = ( + potentialTree + + self.constantTerms(temperature) + ) + + return np.array(potentialTotal) + + def constantTerms(self, temperature: np.ndarray | float) -> np.ndarray | float: + """Need to explicitly compute field-independent but T-dependent parts + that we don't already get from field-dependent loops. At leading order in high-T + expansion these are just (minus) the ideal gas pressure of light particles that + were not integrated over in the one-loop part. + + See Eq. (39) in hep-ph/0510375 for general LO formula + + + Parameters + ---------- + temperature: array-like (float) + The temperature + + Returns + ---------- + constantTerms: array-like (float) + The value of the field-independent contribution to the effective potential + """ + + # How many degrees of freedom we have left. The number of DOFs + # that were included in evaluate() is hardcoded + dofsBoson = self.numBosonDof - 14 + dofsFermion = self.numFermionDof - 12 # we only included top quark loops + + # Fermions contribute with a magic 7/8 prefactor as usual. Overall minus + # sign since Veff(min) = -pressure + return -(dofsBoson + 7.0 / 8.0 * dofsFermion) * np.pi**2 * temperature**4 / 90.0 + + +@pytest.mark.parametrize( + "muSq, lam", + [ + (-50**2, 0.05), + (-50**2, 0.5), + (-100**2, 0.05), + (-100**2, 0.5), + ] +) +def test_EOMSolver(muSq, lam): + """ + Compare the width of the wall computed by WallGo without out-of-equilibrium effects + to the exact result that solves the EOM in the Z2-symmetric quartic potential. + + Parameters + ---------- + muSq : float + Mass term. Must be negative. + lam : float + Quartic coupling. + + Returns + ------- + None. + + """ + + # Exact VEV + vev = np.sqrt(abs(muSq)/lam) + # Exact wall width + trueL = np.sqrt(2/abs(muSq)) + + # Initialise WallGo. Many parameters are now unimportant because we don't solve the + # Boltzmann equation + + Tn = 100 + + temperatureScale = Tn / 100 + fieldScale = vev + derivSettings = WallGo.VeffDerivativeSettings(temperatureVariationScale=float(temperatureScale), fieldValueVariationScale=fieldScale) + + modelParameters = {'muSq': muSq, 'lam': lam} + model = QuarticZ2(modelParameters) + model.getEffectivePotential().configureDerivatives(derivSettings) + + phaseInfo = WallGo.PhaseInfo(temperature = Tn, + phaseLocation1 = WallGo.Fields( [-vev]), + phaseLocation2 = WallGo.Fields( [vev] )) + + thermodynamics = WallGo.Thermodynamics( + model.getEffectivePotential(), + Tn, + phaseInfo.phaseLocation2, + phaseInfo.phaseLocation1, + ) + thermodynamics.freeEnergyHigh.newInterpolationTable(Tn/10, 10*Tn, 100) + thermodynamics.freeEnergyHigh.setExtrapolationType(WallGo.EExtrapolationType.CONSTANT, WallGo.EExtrapolationType.CONSTANT) + thermodynamics.freeEnergyLow.newInterpolationTable(Tn/10, 10*Tn, 100) + thermodynamics.freeEnergyLow.setExtrapolationType(WallGo.EExtrapolationType.CONSTANT, WallGo.EExtrapolationType.CONSTANT) + + # Hack to avoid complains from Hydrodynamics + thermodynamics.TMaxHighT = Tn-1e-6 + thermodynamics.TMaxLowT = Tn-1e-6 + thermodynamics.TMinHighT = Tn-1e-6 + thermodynamics.TMinLowT = Tn-1e-6 + thermodynamics.muMinHighT = 4 + thermodynamics.aMinHighT = 1 + thermodynamics.epsilonMinHighT = 1.5e-3*Tn**4 + thermodynamics.muMaxHighT = 4 + thermodynamics.aMaxHighT = 1 + thermodynamics.epsilonMaxHighT = 1.5e-3*Tn**4 + thermodynamics.muMinLowT = 4 + thermodynamics.aMinLowT = 1-3e-3 + thermodynamics.epsilonMinLowT = 0.0 + thermodynamics.muMaxLowT = 4 + thermodynamics.aMaxLowT = 1-3e-3 + thermodynamics.epsilonMaxLowT = 0.0 + + hydrodynamics = WallGo.Hydrodynamics(thermodynamics, 10, 0.1, rtol=1e-6, atol=1e-8) + grid = WallGo.Grid3Scales( + 50, + 11, + 3*trueL, + 3*trueL, + trueL, + Tn, + 0.75, + 0.1, + ) + boltzmannSolver = WallGo.BoltzmannSolver( + grid, + basisM="Cardinal", + basisN="Chebyshev", + collisionMultiplier=1, + ) + + eom = WallGo.EOM( + boltzmannSolver, + thermodynamics, + hydrodynamics, + grid, + 1, + 3*trueL, + [0.1, 100.0], + [-10,10], + includeOffEq=False, + forceEnergyConservation=False, + forceImproveConvergence=False, + errTol=1e-5, + maxIterations=20, + pressRelErrTol=0.01, + ) + + # Start with a wrong width to see if WallGo can find the correct one. + wallParams = WallGo.WallParams(widths=np.array([1.5*trueL]), offsets=np.array([0])) + # Computes the true wallParams. The wall velocity is irrelevant here. + ( + pressure, + wallParams, + boltzmannResults, + boltzmannBackground, + hydroResults, + *_ + ) = eom.wallPressure(0.3, wallParams) + + tanhError = eom.estimateTanhError(wallParams, boltzmannResults, boltzmannBackground, hydroResults)[0] + + # Compare the exact and WallGo wall widths + np.testing.assert_allclose(trueL, wallParams.widths[0], atol=0, rtol=1e-4) + # Make sure that the tanh error is small + np.testing.assert_allclose(tanhError, 0, atol=1e-8) \ No newline at end of file