From 70f4d7add8f43f658daec67e8279f5acd83352cb Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 5 Feb 2018 12:42:52 +0100
Subject: [PATCH] Bugfixes in LBM phase field code
- 3 phase model works now with step sytem
---
boundaries/boundary_handling.py | 2 +-
lbstep.py | 21 ++--
phasefield/analytical.py | 29 ++---
phasefield/scenario.py | 202 --------------------------------
4 files changed, 23 insertions(+), 231 deletions(-)
delete mode 100644 phasefield/scenario.py
diff --git a/boundaries/boundary_handling.py b/boundaries/boundary_handling.py
index 4157bda3..083f863b 100644
--- a/boundaries/boundary_handling.py
+++ b/boundaries/boundary_handling.py
@@ -32,7 +32,7 @@ class BoundaryHandling:
"If you want to add multiple handlings, choose a different name.")
gpu = self._target == 'gpu'
- dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=gpu)
+ dataHandling.addArray(self._flagFieldName, dtype=np.uint32, cpu=True, gpu=False)
dataHandling.addCustomClass(self._indexArrayName, self.IndexFieldBlockData, cpu=True, gpu=gpu)
ffGhostLayers = self._dataHandling.ghostLayersOfField(self._flagFieldName)
diff --git a/lbstep.py b/lbstep.py
index f07d87f6..66aac26e 100644
--- a/lbstep.py
+++ b/lbstep.py
@@ -42,6 +42,7 @@ class LatticeBoltzmannStep:
Q = len(getStencil(methodParameters['stencil']))
target = optimizationParams['target']
+ self.name = name
self._dataHandling = dataHandling
self._pdfArrName = name + "_pdfSrc"
self._tmpArrName = name + "_pdfTmp"
@@ -67,11 +68,11 @@ class LatticeBoltzmannStep:
if velocityInputArrayName is not None:
methodParameters['velocityInput'] = self._dataHandling.fields[velocityInputArrayName]
- self._kernelParams = kernelParams
+ self.kernelParams = kernelParams
# --- Kernel creation ---
if lbmKernel is None:
- switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, self._kernelParams)
+ switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, self.kernelParams)
optimizationParams['symbolicField'] = dataHandling.fields[self._pdfArrName]
methodParameters['fieldName'] = self._pdfArrName
methodParameters['secondFieldName'] = self._tmpArrName
@@ -82,10 +83,7 @@ class LatticeBoltzmannStep:
self._lbmKernel = lbmKernel
# -- Boundary Handling & Synchronization ---
- if self._gpu:
- self._sync = dataHandling.synchronizationFunctionGPU([self._pdfArrName], methodParameters['stencil'])
- else:
- self._sync = dataHandling.synchronizationFunctionCPU([self._pdfArrName], methodParameters['stencil'])
+ self._sync = dataHandling.synchronizationFunction([self._pdfArrName], methodParameters['stencil'], target)
self._boundaryHandling = BoundaryHandling(self._lbmKernel.method, self._dataHandling, self._pdfArrName,
name=name + "_boundaryHandling",
target=target, openMP=optimizationParams['openMP'])
@@ -98,6 +96,7 @@ class LatticeBoltzmannStep:
for b in self._dataHandling.iterate():
b[self.densityDataName].fill(1.0)
b[self.velocityDataName].fill(0.0)
+ self.setPdfFieldsFromMacroscopicValues()
# -- VTK output
self.vtkWriter = self.dataHandling.vtkWriter(name + str(LatticeBoltzmannStep.vtkScenarioNrCounter),
@@ -180,7 +179,6 @@ class LatticeBoltzmannStep:
return SlicedGetter(self.densitySlice)
def preRun(self):
- self._dataHandling.runKernel(self._setterKernel, **self._kernelParams)
if self._gpu:
self._dataHandling.toGpu(self._pdfArrName)
if self._dataHandling.isOnGpu(self.velocityDataName):
@@ -188,17 +186,20 @@ class LatticeBoltzmannStep:
if self._dataHandling.isOnGpu(self.densityDataName):
self._dataHandling.toGpu(self.densityDataName)
+ def setPdfFieldsFromMacroscopicValues(self):
+ self._dataHandling.runKernel(self._setterKernel, **self.kernelParams)
+
def timeStep(self):
self._sync()
- self._boundaryHandling()
- self._dataHandling.runKernel(self._lbmKernel, **self._kernelParams)
+ self._boundaryHandling(**self.kernelParams)
+ self._dataHandling.runKernel(self._lbmKernel, **self.kernelParams)
self._dataHandling.swap(self._pdfArrName, self._tmpArrName, self._gpu)
self.timeStepsRun += 1
def postRun(self):
if self._gpu:
self._dataHandling.toCpu(self._pdfArrName)
- self._dataHandling.runKernel(self._getterKernel, **self._kernelParams)
+ self._dataHandling.runKernel(self._getterKernel, **self.kernelParams)
def run(self, timeSteps):
self.preRun()
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 1b996d9e..19faba63 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -73,17 +73,8 @@ def symmetricSymbolicSurfaceTension(i, j):
return sp.Symbol("%s_%d_%d" % ((surfaceTensionSymbolName, ) + index))
-def symbolicOrderParameters(numPhases, symbolicLastOrderParameter=False):
- """
- Returns a tuple with numPhases entries, where the all but the last are numbered symbols and the last entry
- is 1 - others
- """
- if symbolicLastOrderParameter:
- return sp.symbols("%s_:%i" % (orderParameterSymbolName, numPhases))
- else:
- phi = sp.symbols("%s_:%i" % (orderParameterSymbolName, numPhases - 1))
- phi = phi + (1 - sum(phi),) # choose last order parameter as 1 - sum(others)
- return phi
+def symbolicOrderParameters(numSymbols):
+ return sp.symbols("%s_:%i" % (orderParameterSymbolName, numSymbols))
def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidthSymbol, transformed=True,
@@ -155,10 +146,12 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
"""
assert not (numPhases is None and orderParameters is None)
if orderParameters is None:
- phi = symbolicOrderParameters(numPhases)
+ phi = symbolicOrderParameters(numPhases-1)
else:
phi = orderParameters
- numPhases = len(phi)
+ numPhases = len(phi) + 1
+
+ phi = tuple(phi) + (1 - sum(phi),)
# Compared to handwritten notes we scale the interface width parameter here to obtain the correct
# equations for the interface profile and the surface tensions i.e. to pass tests
@@ -195,10 +188,10 @@ def analyticInterfaceProfile(x, interfaceWidth=interfaceWidthSymbol):
get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
parameter, i.e. at a transition between two phases.
>>> numPhases = 4
- >>> x, phi = sp.Symbol("x"), symbolicOrderParameters(numPhases)
- >>> F = freeEnergyFunctionalNPhases(numPhases)
+ >>> x, phi = sp.Symbol("x"), symbolicOrderParameters(numPhases-1)
+ >>> F = freeEnergyFunctionalNPhases(orderParameters=phi)
>>> mu = chemicalPotentialsFromFreeEnergy(F)
- >>> mu0 = mu[0].subs({p: 0 for p in phi[1:-1]}) # mu[0] as function of one order parameter only
+ >>> mu0 = mu[0].subs({p: 0 for p in phi[1:]}) # mu[0] as function of one order parameter only
>>> solution = analyticInterfaceProfile(x)
>>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
>>> sp.expand(mu0.subs(solutionSubstitution)) # inserting solution should solve the mu_0=0 equation
@@ -220,7 +213,7 @@ def chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters=None):
return sp.Matrix([functionalDerivative(freeEnergy, op, constants) for op in orderParameters])
-def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx=1):
+def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx=1, cse=True):
"""Reads from order parameter (phi) field and updates chemical potentials"""
chemicalPotential = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
laplaceDiscretization = {Diff(Diff(op)): discreteLaplace(phiField, i, dx)
@@ -229,7 +222,7 @@ def createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiFi
chemicalPotential = chemicalPotential.subs({op: phiField(i) for i, op in enumerate(orderParameters)})
muSweepEqs = [sp.Eq(muField(i), cp) for i, cp in enumerate(chemicalPotential)]
- return muSweepEqs #sympyCseOnEquationList(muSweepEqs)
+ return muSweepEqs if not cse else sympyCseOnEquationList(muSweepEqs)
def createForceUpdateEquations(forceField, phiField, muField, dx=1):
diff --git a/phasefield/scenario.py b/phasefield/scenario.py
deleted file mode 100644
index 8aed3e17..00000000
--- a/phasefield/scenario.py
+++ /dev/null
@@ -1,202 +0,0 @@
-import numpy as np
-
-from lbmpy.boundaries.periodicityhandling import PeriodicityHandling
-from lbmpy.creationfunctions import updateWithDefaultParameters, createLatticeBoltzmannFunction
-from lbmpy.macroscopic_value_kernels import compileMacroscopicValuesSetter
-from lbmpy.phasefield.cahn_hilliard_lbm import createCahnHilliardLbFunction
-from lbmpy.stencils import getStencil
-from lbmpy.updatekernels import createPdfArray
-from pystencils.field import Field, getLayoutOfArray
-from pystencils.slicing import addGhostLayers, removeGhostLayers
-from lbmpy.phasefield.analytical import symbolicOrderParameters, freeEnergyFunctionalNPhases, \
- createChemicalPotentialEvolutionEquations, createForceUpdateEquations
-
-
-class PhasefieldScenario(object):
- def __init__(self, domainSize, numPhases, mobilityRelaxationRates=1.1,
- surfaceTensionCallback=lambda i, j: 1e-3 if i !=j else 0, interfaceWidth=3, dx=1, gamma=1,
- optimizationParams={}, initialVelocity=None, kernelParams={}, **kwargs):
-
- self.numPhases = numPhases
- self.timeStepsRun = 0
- self.domainSize = domainSize
-
- # ---- Parameter normalization
- if not hasattr(mobilityRelaxationRates, '__len__'):
- mobilityRelaxationRates = [mobilityRelaxationRates] * numPhases
-
- D = len(domainSize)
-
- ghostLayers = 1
- domainSizeWithGhostLayer = tuple([s + 2 * ghostLayers for s in domainSize])
-
- if 'stencil' not in kwargs:
- kwargs['stencil'] = 'D2Q9' if D == 2 else 'D3Q27'
-
- methodParameters, optimizationParams = updateWithDefaultParameters(kwargs, optimizationParams)
-
- stencil = getStencil(methodParameters['stencil'])
- fieldLayout = optimizationParams['fieldLayout']
- Q = len(stencil)
-
- if isinstance(initialVelocity, np.ndarray):
- assert initialVelocity.shape[-1] == D
- initialVelocity = addGhostLayers(initialVelocity, indexDimensions=1, ghostLayers=1,
- layout=getLayoutOfArray(self._pdfArrays[0]))
- elif initialVelocity is None:
- initialVelocity = [0] * D
-
- self.kernelParams = kernelParams
-
- # ---- Arrays
- self.velArr = np.zeros(domainSizeWithGhostLayer + (D,), order=fieldLayout)
- self.muArr = np.zeros(domainSizeWithGhostLayer + (numPhases - 1,), order=fieldLayout)
- self.phiArr = np.zeros(domainSizeWithGhostLayer + (numPhases - 1,), order=fieldLayout)
- self.forceArr = np.zeros(domainSizeWithGhostLayer + (D,), order=fieldLayout)
-
- self._pdfArrays = [[createPdfArray(domainSize, Q, layout=optimizationParams['fieldLayout'])
- for i in range(numPhases)],
- [createPdfArray(domainSize, Q, layout=optimizationParams['fieldLayout'])
- for i in range(numPhases)]]
-
- # ---- Fields
- velField = Field.createFromNumpyArray('vel', self.velArr, indexDimensions=1)
- muField = Field.createFromNumpyArray('mu', self.muArr, indexDimensions=1)
- phiField = Field.createFromNumpyArray('phi', self.phiArr, indexDimensions=1)
- forceField = Field.createFromNumpyArray('F', self.forceArr, indexDimensions=1)
-
- orderParameters = symbolicOrderParameters(numPhases)
- freeEnergy = freeEnergyFunctionalNPhases(numPhases, surfaceTensionCallback, interfaceWidth, orderParameters)
-
- # ---- Sweeps
- muSweepEquations = createChemicalPotentialEvolutionEquations(freeEnergy, orderParameters, phiField, muField, dx)
- forceSweepEquations = createForceUpdateEquations(numPhases, forceField, phiField, muField, dx)
- if optimizationParams['target'] == 'cpu':
- from pystencils.cpu import createKernel, makePythonFunction
- self.muSweep = makePythonFunction(createKernel(muSweepEquations))
- self.forceSweep = makePythonFunction(createKernel(forceSweepEquations))
- else:
- from pystencils.gpucuda import createCUDAKernel, makePythonFunction
- self.muSweep = makePythonFunction(createCUDAKernel(muSweepEquations))
- self.forceSweep = makePythonFunction(createCUDAKernel(forceSweepEquations))
-
- optimizationParams['pdfArr'] = self._pdfArrays[0][0]
-
- self.lbSweepHydro = createLatticeBoltzmannFunction(force=[forceField(i) for i in range(D)],
- output={'velocity': velField},
- optimizationParams=optimizationParams, **kwargs)
-
- useFdForCahnHilliard = False
- 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'])
-
- assert self.muArr.shape == self.phiArr.shape
- self._muPhiPeriodicityHandler = PeriodicityHandling(self.muArr.shape, (True, True, True),
- optimizationParams['target'])
-
- # Pdf array initialization
- hydroLbmInit = compileMacroscopicValuesSetter(self.lbSweepHydro.method,
- {'density': 1.0, 'velocity': initialVelocity},
- pdfArr=self._pdfArrays[0][0], target='cpu')
- hydroLbmInit(pdfs=self._pdfArrays[0][0], F=self.forceArr, **self.kernelParams)
- self.initializeCahnHilliardPdfsAccordingToPhi()
-
- self._nonPdfArrays = {
- 'phiArr': self.phiArr,
- 'muArr': self.muArr,
- 'velArr': self.velArr,
- 'forceArr': self.forceArr,
- }
- self._nonPdfGpuArrays = None
- self._pdfGpuArrays = None
- self.target = optimizationParams['target']
-
- self.hydroVelocitySetter = None
-
- def updateHydroPdfsAccordingToVelocity(self):
- if self.hydroVelocitySetter is None:
- self.hydroVelocitySetter = compileMacroscopicValuesSetter(self.lbSweepHydro.method,
- {'density': 1.0, 'velocity': self.velArr},
- pdfArr=self._pdfArrays[0][0], target='cpu')
- self.hydroVelocitySetter(pdfs=self._pdfArrays[0][0], F=self.forceArr, **self.kernelParams)
-
- def _arraysFromCpuToGpu(self):
- import pycuda.gpuarray as gpuarray
- if self._nonPdfGpuArrays is None:
- self._nonPdfGpuArrays = {name: gpuarray.to_gpu(arr) for name, arr in self._nonPdfArrays.items()}
- self._pdfGpuArrays = [[gpuarray.to_gpu(arr) for arr in self._pdfArrays[0]],
- [gpuarray.to_gpu(arr) for arr in self._pdfArrays[1]]]
- else:
- for name, arr in self._nonPdfArrays.items():
- self._nonPdfGpuArrays[name].set(arr)
- for i in range(2):
- for cpuArr, gpuArr in zip(self._pdfArrays[i], self._pdfGpuArrays[i]):
- gpuArr.set(cpuArr)
-
- def _arraysFromGpuToCpu(self):
- for name, arr in self._nonPdfArrays.items():
- self._nonPdfGpuArrays[name].get(arr)
- for cpuArr, gpuArr in zip(self._pdfArrays[0], self._pdfGpuArrays[0]):
- gpuArr.get(cpuArr)
-
- def initializeCahnHilliardPdfsAccordingToPhi(self):
- for i in range(1, self.numPhases):
- self._pdfArrays[0][i].fill(0)
- self._pdfArrays[0][i][..., 0] = self.phiArr[..., i-1]
-
- def gaussianSmoothPhiFields(self, sigma):
- from scipy.ndimage.filters import gaussian_filter
- for i in range(self.phiArr.shape[-1]):
- gaussian_filter(self.phi[..., i], sigma, output=self.phi[..., i], mode='wrap')
-
- @property
- def phi(self):
- return removeGhostLayers(self.phiArr, indexDimensions=1)
-
- @property
- def mu(self):
- return removeGhostLayers(self.muArr, indexDimensions=1)
-
- @property
- def velocity(self):
- return removeGhostLayers(self.velArr, indexDimensions=1)
-
- def run(self, timeSteps=1):
- """Run the scenario for the given amount of time steps"""
- if self.target == 'gpu':
- self._arraysFromCpuToGpu()
- self._timeLoop(self._pdfGpuArrays, timeSteps=timeSteps, **self._nonPdfGpuArrays)
- self._arraysFromGpuToCpu()
- else:
- self._timeLoop(self._pdfArrays, timeSteps=timeSteps, **self._nonPdfArrays)
- self.timeStepsRun += timeSteps
-
- def _timeLoop(self, pdfArrays, phiArr, muArr, velArr, forceArr, timeSteps):
- for t in range(timeSteps):
- self._muPhiPeriodicityHandler(pdfs=phiArr)
- self.muSweep(phi=phiArr, mu=muArr)
-
- self._muPhiPeriodicityHandler(pdfs=muArr)
- self.forceSweep(mu=muArr, phi=phiArr, F=forceArr)
-
- for src in pdfArrays[0]:
- self._pdfPeriodicityHandler(pdfs=src)
-
- for sweep, src, dst in zip(self.lbSweeps, *pdfArrays):
- sweep(src=src, dst=dst, F=forceArr, phi=phiArr, vel=velArr, mu=muArr)
-
- pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0]
--
GitLab