diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index ca84f6b53b14ede4f96a3779c52ba906150c44c5..d98ac5423e859498749c42507ba086acfedf1ec4 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -29,8 +29,8 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt
kappaPrime = tuple(interfaceWidth**2 * k for k in kappa)
C = sp.symbols("C_:3")
- bulkFreeEnergy = sum(k / 2 * C_i ** 2 * (1 - C_i) ** 2 for k, C_i in zip(kappa, C))
- surfaceFreeEnergy = sum(k / 2 * Diff(C_i) ** 2 for k, C_i in zip(kappaPrime, C))
+ bulkFreeEnergy = sum(k * C_i ** 2 * (1 - C_i) ** 2 / 2 for k, C_i in zip(kappa, C))
+ surfaceFreeEnergy = sum(k * Diff(C_i) ** 2 / 2 for k, C_i in zip(kappaPrime, C))
F = 0
if includeBulk:
@@ -59,6 +59,25 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt
return F
+def freeEnergyFunctionalNPhasesPenaltyTerm(orderParameters, interfaceWidth=interfaceWidthSymbol, kappa=None,
+ penaltyTermFactor=0.01):
+ numPhases = len(orderParameters)
+ if kappa is None:
+ kappa = sp.symbols("kappa_:%d" % (numPhases,))
+ if not hasattr(kappa, "__length__"):
+ kappa = [kappa] * numPhases
+
+ def f(x):
+ return x ** 2 * (1 - x) ** 2
+
+ bulk = sum(f(c) * k / 2 for c, k in zip(orderParameters, kappa))
+ interface = sum(Diff(c) ** 2 / 2 * interfaceWidth ** 2 * k
+ for c, k in zip(orderParameters, kappa))
+
+ bulkPenaltyTerm = (1 - sum(c for c in orderParameters)) ** 2
+ return bulk + interface + penaltyTermFactor * bulkPenaltyTerm
+
+
def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymbolicSurfaceTension,
interfaceWidth=interfaceWidthSymbol, orderParameters=None,
includeBulk=True, includeInterface=True, symbolicLambda=False,
diff --git a/phasefield/experiments1D.py b/phasefield/experiments1D.py
index ae304d6c0f9b9c494aec276f3eea71ec99e1058f..867e09024034113efec2097e6171131a40a5ea02 100644
--- a/phasefield/experiments1D.py
+++ b/phasefield/experiments1D.py
@@ -1,4 +1,7 @@
import numpy as np
+import sympy as sp
+
+from lbmpy.chapman_enskog import Diff
from pystencils import makeSlice
@@ -33,24 +36,132 @@ def plotStatus(phaseFieldStep, fromX=None, toX=None):
plt.legend()
-def initSharpInterface(pfStep, phaseIdx1=1, phaseIdx2=2, x1=None, x2=None):
+def plotFreeEnergyBulkContours(freeEnergy, orderParameters, phase0=0, phase1=1, **kwargs):
+ import lbmpy.plot2d as plt
+
+ x = np.linspace(-.2, 1.2, 100)
+ y = np.linspace(-.2, 1.2, 100)
+ xg, yg = np.meshgrid(x, y)
+ substitutions = {op: 0 for i, op in enumerate(orderParameters) if i not in (phase0, phase1)}
+ substitutions.update({d: 0 for d in freeEnergy.atoms(Diff)}) # remove interface components of free energy
+ freeEnergyLambda = sp.lambdify([orderParameters[phase0], orderParameters[phase1]],
+ freeEnergy.subs(substitutions))
+ if 'levels' not in kwargs:
+ kwargs['levels'] = np.linspace(0, 1, 60)
+ plt.contour(x, y, freeEnergyLambda(xg, yg), **kwargs)
+
+
+def initSharpInterface(pfStep, phaseIdx, x1=None, x2=None, inverse=False):
domainSize = pfStep.dataHandling.shape
if x1 is None:
x1 = domainSize[0] // 4
if x2 is None:
x2 = 3 * x1
+ if phaseIdx >= pfStep.numOrderParameters:
+ return
+
for b in pfStep.dataHandling.iterate():
x = b.cellIndexArrays[0]
mid = np.logical_and(x1 < x, x < x2)
phi = b[pfStep.phiFieldName]
+ val1, val2 = (1, 0) if inverse else (0, 1)
+
+ phi[..., phaseIdx].fill(val1)
+ phi[mid, phaseIdx] = val2
+
+ pfStep.setPdfFieldsFromMacroscopicValues()
- if phaseIdx1 is not None:
- phi[..., phaseIdx1].fill(0)
- phi[mid, phaseIdx1] = 1
- if phaseIdx2 is not None:
- phi[..., phaseIdx2].fill(1)
- phi[mid, phaseIdx2] = 0
+def tanhTest(pfStep, phase0, phase1, expectedInterfaceWidth=1, timeSteps=10000):
+ """
+ Initializes a sharp interface and checks if tanh-shaped profile is developing
+ :param pfStep: phase field scenario / step
+ :param phase0: index of first phase to initialize
+ :param phase1: index of second phase to initialize inversely
+ :param expectedInterfaceWidth: interface width parameter alpha that is used in analytical form
+ :param timeSteps: number of time steps run before evaluation
+ :return: deviation of simulated profile from analytical solution as average(abs(simulation-analytic))
+ """
+ import lbmpy.plot2d as plt
+ from lbmpy.phasefield.analytical import analyticInterfaceProfile
+
+ domainSize = pfStep.dataHandling.shape
+ pfStep.reset()
+ pfStep.dataHandling.fill(pfStep.phiFieldName, 0)
+ initSharpInterface(pfStep, phaseIdx=phase0, inverse=False)
+ initSharpInterface(pfStep, phaseIdx=phase1, inverse=True)
pfStep.setPdfFieldsFromMacroscopicValues()
+ pfStep.run(timeSteps)
+
+ visWidth = 20
+ x = np.arange(visWidth) - (visWidth // 2)
+ analytic = np.array([analyticInterfaceProfile(x_i - 0.5, expectedInterfaceWidth) for x_i in x], dtype=np.float64)
+
+ stepLocation = domainSize[0] // 4
+ simulated = pfStep.phi[stepLocation - visWidth//2:stepLocation + visWidth//2, 0, phase0]
+ plt.plot(analytic, label='analytic', marker='o')
+ plt.plot(simulated, label='simulated', marker='x')
+ plt.legend()
+
+ return np.average(np.abs(simulated - analytic))
+
+
+def galileanInvarianceTest(pfStep, velocity=0.05, rounds=3, phase0=0, phase1=1,
+ expectedInterfaceWidth=1, initTimeSteps=5000):
+ """
+ Moves interface at constant speed through periodic domain - check if momentum is conserved
+ :param pfStep: phase field scenario / step
+ :param velocity: constant velocity to move interface
+ :param rounds: how many times the interface should travel through the domain
+ :param phase0: index of first phase to initialize
+ :param phase1: index of second phase to initialize inversely
+ :param expectedInterfaceWidth: interface width parameter alpha that is used in analytical form
+ :param initTimeSteps: before velocity is set, this many time steps are run to let interface settle to tanh shape
+ :return: change in velocity
+ """
+ import lbmpy.plot2d as plt
+ from lbmpy.phasefield.analytical import analyticInterfaceProfile
+
+ domainSize = pfStep.dataHandling.shape
+ roundTimeSteps = int((domainSize[0]+0.25) / velocity)
+
+ print("Velocity:", velocity, " Timesteps for round:", roundTimeSteps)
+
+ pfStep.reset()
+ pfStep.dataHandling.fill(pfStep.phiFieldName, 0)
+ initSharpInterface(pfStep, phaseIdx=phase0, inverse=False)
+ initSharpInterface(pfStep, phaseIdx=phase1, inverse=True)
+ pfStep.setPdfFieldsFromMacroscopicValues()
+
+ print("Running", initTimeSteps, "initial time steps")
+ pfStep.run(initTimeSteps)
+ pfStep.dataHandling.fill(pfStep.velFieldName, velocity, fValue=0)
+ pfStep.setPdfFieldsFromMacroscopicValues()
+
+ stepLocation = domainSize[0] // 4
+ visWidth = 20
+
+ simulatedProfiles = []
+
+ def captureProfile():
+ simulated = pfStep.phi[stepLocation - visWidth // 2:stepLocation + visWidth // 2, 0, phase0].copy()
+ simulatedProfiles.append(simulated)
+
+ captureProfile()
+ for rt in range(rounds ):
+ print("Running round %d/%d" % (rt+1, rounds))
+ pfStep.run(roundTimeSteps)
+ captureProfile()
+
+ x = np.arange(visWidth) - (visWidth // 2)
+ analytic = np.array([analyticInterfaceProfile(x_i - 0.5, expectedInterfaceWidth) for x_i in x], dtype=np.float64)
+
+ plt.plot(x, analytic, label='analytic', marker='o')
+ for i, profile in enumerate(simulatedProfiles):
+ plt.plot(x, profile, label="After %d rounds" % (i,))
+
+ plt.legend()
+
+ return np.average(pfStep.velocity[:, 0, 0]) - velocity
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index e6dc213dcbddd75d038cfae7e152ad66027359cb..98245be4b68c91aea808d23204ea18f05d9530e3 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -17,7 +17,7 @@ class PhaseFieldStep:
target='cpu', openMP=False, kernelParams={}, dx=1, dt=1, solveCahnHilliardWithFiniteDifferences=False):
if dataHandling is None:
- dataHandling = SerialDataHandling(domainSize)
+ dataHandling = SerialDataHandling(domainSize, periodicity=True)
self.freeEnergy = freeEnergy
self.chemicalPotentials = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
@@ -95,17 +95,22 @@ class PhaseFieldStep:
name=name + "_chLbm_%d" % (i,), )
self.cahnHilliardSteps.append(chStep)
+ self.runHydroLbm = True
+ self.densityOrderParameter = densityOrderParameter
+ self.timeStepsRun = 0
+ self.reset()
+
+ def reset(self):
# Init φ and μ
self.dataHandling.fill(self.phiFieldName, 0.0)
- self.dataHandling.fill(self.phiFieldName, 1.0 if densityOrderParameter is not None else 0.0, fValue=0)
+ self.dataHandling.fill(self.phiFieldName, 1.0 if self.densityOrderParameter is not None else 0.0, fValue=0)
self.dataHandling.fill(self.muFieldName, 0.0)
self.dataHandling.fill(self.forceFieldName, 0.0)
+ self.dataHandling.fill(self.velFieldName, 0.0)
self.setPdfFieldsFromMacroscopicValues()
self.timeStepsRun = 0
- self.runHydroLbm = True
-
def setPdfFieldsFromMacroscopicValues(self):
self.hydroLbmStep.setPdfFieldsFromMacroscopicValues()
for chStep in self.cahnHilliardSteps: