diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index a2e0f5f7a34ed22899c5377daaa517133d441701..5cae60b6ea8b261eacc9fc6fb6b7ee8a2d4a6b55 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -177,3 +177,17 @@ def createForceUpdateEquations(numPhases, forceField, phiField, muField, dx=1):
forceSweepEqs.append(sp.Eq(forceField(d), rhs))
return forceSweepEqs
+
+def cahnHilliardFdKernel(phaseIdx, phi, mu, velocity, mobility, dx, dt, target='cpu'):
+ from pystencils.finitedifferences import transient, advection, diffusion, Discretization2ndOrder
+ cahnHilliard = transient(phi, phaseIdx) + advection(phi, velocity, phaseIdx) - diffusion(mu, mobility, phaseIdx)
+ discretizedEq = Discretization2ndOrder(dx, dt)(cahnHilliard)
+ updateRule = [sp.Eq(phi.newFieldWithDifferentName('dst')(phaseIdx),
+ discretizedEq)]
+
+ if target == 'cpu':
+ from pystencils.cpu import createKernel, makePythonFunction
+ return makePythonFunction(createKernel(updateRule))
+ elif target == 'gpu':
+ from pystencils.gpucuda import createCUDAKernel, makePythonFunction
+ return makePythonFunction(createCUDAKernel(updateRule))
diff --git a/phasefield/plot.py b/phasefield/plot.py
index eb58f368fb9413ebd7366aa257cf7f8c9c1c0c8b..652943f9bb3435745a01e7a33951c3e625ee61fa 100644
--- a/phasefield/plot.py
+++ b/phasefield/plot.py
@@ -1,7 +1,7 @@
from lbmpy.plot2d import *
-def drawAngles(ax, groupedPoints):
+def _drawAngles(ax, groupedPoints):
from matplotlib.lines import Line2D
from matplotlib.text import Annotation
from lbmpy.phasefield.postprocessing import getAngle
@@ -32,25 +32,33 @@ def drawAngles(ax, groupedPoints):
ax.add_line(lines)
angle = getAngle(*groupedPoints)
- text = Annotation('%.1f°' % (angle,), (textPos[1], textPos[0]), axes=ax, fontsize=15, ha=ha, va=va)
- ax.add_artist(text)
+ annotation = Annotation('%.1f°' % (angle,), (textPos[1], textPos[0]), axes=ax, fontsize=15, ha=ha, va=va)
+ ax.add_artist(annotation)
-def plotAngles(phi0, phi1, phi2, branchingDistance=0.5, onlyFirst=True):
+def phaseIndices(phi, **kwargs):
+ singlePhaseArrays = [phi[..., i] for i in range(phi.shape[-1])]
+ lastPhase = 1 - sum(singlePhaseArrays)
+ singlePhaseArrays.append(lastPhase)
+ idxArray = np.array(singlePhaseArrays).argmax(axis=0)
+ return scalarField(idxArray, **kwargs)
+
+
+def angles(phi0, phi1, phi2, branchingDistance=0.5, branchingPointFilter=3, onlyFirst=True):
from lbmpy.phasefield.postprocessing import getTriplePointInfos
levels = [0.5, ]
scalarFieldContour(phi0, levels=levels)
scalarFieldContour(phi1, levels=levels)
- scalarFieldContour(phi2, levels=levels);
+ scalarFieldContour(phi2, levels=levels)
if onlyFirst:
- angleInfo = getTriplePointInfos(phi0, phi1, phi2, branchingDistance)
+ angleInfo = getTriplePointInfos(phi0, phi1, phi2, branchingDistance, branchingPointFilter)
else:
- angleInfo = getTriplePointInfos(phi0, phi1, phi2, branchingDistance) +\
- getTriplePointInfos(phi1, phi0, phi2, branchingDistance) +\
- getTriplePointInfos(phi2, phi1, phi0, branchingDistance)
+ angleInfo = getTriplePointInfos(phi0, phi1, phi2, branchingDistance, branchingPointFilter) +\
+ getTriplePointInfos(phi1, phi0, phi2, branchingDistance, branchingPointFilter) +\
+ getTriplePointInfos(phi2, phi1, phi0, branchingDistance, branchingPointFilter)
for points in angleInfo:
- drawAngles(gca(), points)
+ _drawAngles(gca(), points)
diff --git a/phasefield/postprocessing.py b/phasefield/postprocessing.py
index 63debd11196be656fb7f677348c1df676c3d8a5b..a6c63bea0f652bf8867435166a75fd6e5d1b6f72 100644
--- a/phasefield/postprocessing.py
+++ b/phasefield/postprocessing.py
@@ -33,19 +33,19 @@ def findJumpIndices(array, threshold=0, minLength=3):
array = array[jump + minLength:]
-def findBranchingPoint(pathVertices1, pathVertices2, maxDistance=0.5):
+def findBranchingPoint(pathVertices1, pathVertices2, maxDistance=0.5, branchingPointFilter=3):
tree = scipy.spatial.KDTree(pathVertices1)
distances, indices = tree.query(pathVertices2, k=1, distance_upper_bound=maxDistance)
distances[distances == np.inf] = -1
- jumpIndices = findJumpIndices(distances, 0, 3)
+ jumpIndices = findJumpIndices(distances, 0, branchingPointFilter)
return pathVertices2[jumpIndices]
-def findAllBranchingPoints(phaseField1, phaseField2, maxDistance=0.1):
+def findAllBranchingPoints(phaseField1, phaseField2, maxDistance=0.1, branchingPointFilter=3):
result = []
isoLines = [getIsolines(p, level=0.5, refinementFactor=4) for p in (phaseField1, phaseField2)]
for path1, path2 in itertools.product(*isoLines):
- bbs = findBranchingPoint(path1, path2, maxDistance)
+ bbs = findBranchingPoint(path1, path2, maxDistance, branchingPointFilter)
result += list(bbs)
return np.array(result)
@@ -114,7 +114,7 @@ def getAngle(pMid, p1, p2):
return result
-def getTriplePointInfos(phi0, phi1, phi2, branchingDistance=0.5):
+def getTriplePointInfos(phi0, phi1, phi2, branchingDistance=0.5, branchingPointFilter=3):
"""
:param branchingDistance: where the 1/2 contour lines move apart farther than this value, the
@@ -126,9 +126,11 @@ def getTriplePointInfos(phi0, phi1, phi2, branchingDistance=0.5):
# the angle at the triple points is measured with contour lines of level 1/2 at "branching points"
# i.e. at points where the lines move away from each other
- bb1 = findAllBranchingPoints(phi0, phi1, branchingDistance)
- bb2 = findAllBranchingPoints(phi0, phi2, branchingDistance)
+ bb1 = findAllBranchingPoints(phi0, phi1, branchingDistance, branchingPointFilter)
+ bb2 = findAllBranchingPoints(phi0, phi2, branchingDistance, branchingPointFilter)
ip = findAllIntersectionPoints(phi0, phi1)
return groupPoints(ip, np.vstack([bb1, bb2]))
+def getAnglesOfPhase(phi0, phi1, phi2, branchingDistance=0.5, branchingPointFilter=3):
+ return [getAngle(*r) for r in getTriplePointInfos(phi0, phi1, phi2, branchingDistance, branchingPointFilter)]
diff --git a/phasefield/scenario.py b/phasefield/scenario.py
index 2d55802d735fd5af40891cc2bf42dad5467e2a9e..97eed47f2cd3591e06abffefb5511a049ea656a9 100644
--- a/phasefield/scenario.py
+++ b/phasefield/scenario.py
@@ -85,11 +85,21 @@ class PhasefieldScenario(object):
self.lbSweepHydro = createLatticeBoltzmannFunction(force=[forceField(i) for i in range(D)],
output={'velocity': velField},
optimizationParams=optimizationParams, **kwargs)
- self.lbSweepsCH = [createCahnHilliardLbFunction(stencil, mobilityRelaxationRates[i],
- velField, muField(i), phiField(i), optimizationParams, gamma)
- for i in range(numPhases-1)]
- self.lbSweeps = [self.lbSweepHydro] + self.lbSweepsCH
+ useFdForCahnHilliard = True
+ if useFdForCahnHilliard:
+ dt = 0.01
+ mobility = 1
+ from lbmpy.phasefield.analytical import cahnHilliardFdKernel
+ self.sweepsCH = [cahnHilliardFdKernel(i, phiField, muField, velField, mobility,
+ dx, dt, optimizationParams['target'])
+ for i in range(numPhases-1)]
+ else:
+ self.sweepsCH = [createCahnHilliardLbFunction(stencil, mobilityRelaxationRates[i],
+ velField, muField(i), phiField(i), optimizationParams, gamma)
+ for i in range(numPhases-1)]
+
+ self.lbSweeps = [self.lbSweepHydro] + self.sweepsCH
self._pdfPeriodicityHandler = PeriodicityHandling(self._pdfArrays[0][0].shape, (True, True, True),
optimizationParams['target'])