Commit 9cf1ac28 authored by Martin Bauer's avatar Martin Bauer
Browse files

lbmpy phasefield

- step class for LB phasefield generic enough to work with 3-phase and
  N-phase models
- cahn hilliard can either be solved by LBM or by finite differences
- 3 phase model can be solved with rho phase or without
parent 4e7953f5
import sympy as sp
from pystencils.bitoperations import xor, rightShift, leftShift
try:
from sympy.utilities.codegen import CCodePrinter
except ImportError:
......@@ -210,6 +213,12 @@ class CustomSympyPrinter(CCodePrinter):
if expr.func == castFunc:
arg, type = expr.args
return "*((%s)(& %s))" % (PointerType(type), self._print(arg))
elif expr.func == xor:
return "(%s ^ %s" % (self._print(expr.args[0]), self._print(expr.args[1]))
elif expr.func == rightShift:
return "(%s >> %s)" % (self._print(expr.args[0]), self._print(expr.args[1]))
elif expr.func == leftShift:
return "(%s << %s)" % (self._print(expr.args[0]), self._print(expr.args[1]))
else:
return super(CustomSympyPrinter, self)._print_Function(expr)
......
import sympy as sp
xor = sp.Function("⊻")
rightShift = sp.Function("rshift")
leftShift = sp.Function("lshift")
......@@ -241,12 +241,13 @@ class DataHandling(ABC):
def __str__(self):
result = ""
rowFormat = "{:>35}|{:>21}|{:>21}\n"
separatorLine = "" * (25 + 21 + 21 + 2) + "\n"
firstColumnWidth = max(len("Name"), max(len(a) for a in self.arrayNames))
rowFormat = "{:>%d}|{:>21}|{:>21}\n" % (firstColumnWidth,)
separatorLine = "-" * (firstColumnWidth + 21 + 21 + 2) + "\n"
result += rowFormat.format("Name", "Inner (min/max)", "WithGl (min/max)")
result += separatorLine
for arrName in sorted(self.arrayNames):
innerMinMax = (self.min(arrName, ghostLayers=False), self.max(arrName, ghostLayers=False))
innerMinMax = (self.min(arrName, ghostLayers=False), self.max(arrName, ghostLayers=False))
withGlMinMax = (self.min(arrName, ghostLayers=True), self.max(arrName, ghostLayers=True))
innerMinMax = "({0[0]:3.3g},{0[1]:3.3g})".format(innerMinMax)
withGlMinMax = "({0[0]:3.3g},{0[1]:3.3g})".format(withGlMinMax)
......
......@@ -84,8 +84,6 @@ class ParallelDataHandling(DataHandling):
ghostLayers = self.defaultGhostLayers
if layout is None:
layout = self.defaultLayout
if latexName is None:
latexName = name
if len(self.blocks) == 0:
raise ValueError("Data handling expects that each process has at least one block")
if hasattr(dtype, 'type'):
......@@ -118,13 +116,14 @@ class ParallelDataHandling(DataHandling):
if indexDimensions == 1:
shape += (fSize, )
assert all(f.name != latexName for f in self.fields.values()), "Symbolic field with this name already exists"
assert all(f.name != name for f in self.fields.values()), "Symbolic field with this name already exists"
self.fields[name] = Field.createGeneric(latexName, self.dim, dtype, indexDimensions, layout,
self.fields[name] = Field.createGeneric(name, self.dim, dtype, indexDimensions, layout,
indexShape=(fSize,) if indexDimensions > 0 else None)
self._fieldNameToCpuDataName[latexName] = name
self.fields[name].latexName = latexName
self._fieldNameToCpuDataName[name] = name
if gpu:
self._fieldNameToGpuDataName[latexName] = self.GPU_DATA_PREFIX + name
self._fieldNameToGpuDataName[name] = self.GPU_DATA_PREFIX + name
return self.fields[name]
def hasData(self, name):
......@@ -139,7 +138,7 @@ class ParallelDataHandling(DataHandling):
return tuple(self._customDataNames)
def addArrayLike(self, name, nameOfTemplateField, latexName=None, cpu=True, gpu=False):
self.addArray(name, latexName=latexName, cpu=cpu, gpu=gpu, **self._fieldInformation[nameOfTemplateField])
return self.addArray(name, latexName=latexName, cpu=cpu, gpu=gpu, **self._fieldInformation[nameOfTemplateField])
def swap(self, name1, name2, gpu=False):
if gpu:
......@@ -178,16 +177,21 @@ class ParallelDataHandling(DataHandling):
if sliceObj is None:
sliceObj = tuple([slice(None, None, None)] * self.dim)
if self.dim == 2:
sliceObj += (0.5,)
sliceObj = sliceObj[:2] + (0.5,) + sliceObj[2:]
array = wlb.field.gatherField(self.blocks, name, sliceObj, allGather)
lastElement = sliceObj[3:]
array = wlb.field.gatherField(self.blocks, name, sliceObj[:3], allGather)
if array is None:
return None
if self.fields[name].indexDimensions == 0:
array = array[..., 0]
if self.dim == 2:
array = array[:, :, 0]
if lastElement and self.fields[name].indexDimensions > 0:
array = array[..., lastElement[0]]
if self.fields[name].indexDimensions == 0:
array = array[..., 0]
return array
def _normalizeArrShape(self, arr, indexDimensions):
......
......@@ -31,7 +31,6 @@ class SerialDataHandling(DataHandling):
self.defaultGhostLayers = defaultGhostLayers
self.defaultLayout = defaultLayout
self._fields = DotDict()
self._fieldLatexNameToDataName = {}
self.cpuArrays = DotDict()
self.gpuArrays = DotDict()
self.customDataCpu = DotDict()
......@@ -74,8 +73,6 @@ class SerialDataHandling(DataHandling):
ghostLayers = self.defaultGhostLayers
if layout is None:
layout = self.defaultLayout
if latexName is None:
latexName = name
kwargs = {
'shape': tuple(s + 2 * ghostLayers for s in self._domainSize),
......@@ -109,10 +106,10 @@ class SerialDataHandling(DataHandling):
raise ValueError("GPU Field with this name already exists")
self.gpuArrays[name] = gpuarray.to_gpu(cpuArr)
assert all(f.name != latexName for f in self.fields.values()), "Symbolic field with this name already exists"
self.fields[name] = Field.createFixedSize(latexName, shape=kwargs['shape'], indexDimensions=indexDimensions,
assert all(f.name != name for f in self.fields.values()), "Symbolic field with this name already exists"
self.fields[name] = Field.createFixedSize(name, shape=kwargs['shape'], indexDimensions=indexDimensions,
dtype=kwargs['dtype'], layout=layoutTuple)
self._fieldLatexNameToDataName[latexName] = name
self.fields[name].latexName = latexName
return self.fields[name]
def addCustomData(self, name, cpuCreationFunction,
......@@ -136,7 +133,7 @@ class SerialDataHandling(DataHandling):
return name in self.fields
def addArrayLike(self, name, nameOfTemplateField, latexName=None, cpu=True, gpu=False):
self.addArray(name, latexName=latexName, cpu=cpu, gpu=gpu, **self._fieldInformation[nameOfTemplateField])
return self.addArray(name, latexName=latexName, cpu=cpu, gpu=gpu, **self._fieldInformation[nameOfTemplateField])
def iterate(self, sliceObj=None, gpu=False, ghostLayers=True, innerGhostLayers=True):
if ghostLayers is True:
......@@ -172,12 +169,15 @@ class SerialDataHandling(DataHandling):
glToRemove = 0
arr = self.cpuArrays[name]
indDimensions = self.fields[name].indexDimensions
spatialDimensions = self.fields[name].spatialDimensions
arr = removeGhostLayers(arr, indexDimensions=indDimensions, ghostLayers=glToRemove)
if sliceObj is not None:
sliceObj = normalizeSlice(sliceObj, arr.shape[:-indDimensions] if indDimensions > 0 else arr.shape)
sliceObj = tuple(s if type(s) is slice else slice(s, s + 1, None) for s in sliceObj)
arr = arr[sliceObj]
normalizedSlice = normalizeSlice(sliceObj[:spatialDimensions], arr.shape[:spatialDimensions])
normalizedSlice = tuple(s if type(s) is slice else slice(s, s + 1, None) for s in normalizedSlice)
normalizedSlice += sliceObj[spatialDimensions:]
arr = arr[normalizedSlice]
else:
arr = arr.view()
arr.flags.writeable = False
......@@ -198,7 +198,7 @@ class SerialDataHandling(DataHandling):
self.toGpu(name)
def runKernel(self, kernelFunc, *args, **kwargs):
dataUsedInKernel = [self._fieldLatexNameToDataName[p.fieldName]
dataUsedInKernel = [p.fieldName
for p in kernelFunc.parameters if p.isFieldPtrArgument and p.fieldName not in kwargs]
arrays = self.gpuArrays if kernelFunc.ast.backend == 'gpucuda' else self.cpuArrays
arrayParams = {name: arrays[name] for name in dataUsedInKernel}
......@@ -270,11 +270,11 @@ class SerialDataHandling(DataHandling):
if target == 'cpu':
def resultFunctor():
for func in resultFunctors:
for name, func in zip(names, resultFunctors):
func(pdfs=self.cpuArrays[name])
else:
def resultFunctor():
for func in resultFunctors:
for name, func in zip(names, resultFunctors):
func(pdfs=self.gpuArrays[name])
return resultFunctor
......
......@@ -188,6 +188,7 @@ class Field(object):
self._layout = normalizeLayout(layout)
self.shape = shape
self.strides = strides
self.latexName = None
def newFieldWithDifferentName(self, newName):
return Field(newName, self.fieldType, self._dtype, self._layout, self.shape, self.strides)
......@@ -381,10 +382,11 @@ class Field(object):
return self._offsetName
def _latex(self, arg):
n = self._field.latexName if self._field.latexName else self._field.name
if self._superscript:
return "{{%s}_{%s}^{%s}}" % (self._field.name, self._offsetName, self._superscript)
return "{{%s}_{%s}^{%s}}" % (n, self._offsetName, self._superscript)
else:
return "{{%s}_{%s}}" % (self._field.name, self._offsetName)
return "{{%s}_{%s}}" % (n, self._offsetName)
@property
def index(self):
......
......@@ -247,7 +247,7 @@ class Transient(sp.Function):
@property
def scalar(self):
if self.scalarIndex is None:
return self.args[0].field
return self.args[0].field.center
else:
return self.args[0].field(self.scalarIndex)
......@@ -310,9 +310,9 @@ class Discretization2ndOrder:
if len(solveResult) != 1:
raise ValueError("Could not solve for transient term" + str(solveResult))
rhs = solveResult.pop()
idx = 0 if transientTerm.scalarIndex is None else transientTerm.scalarIndex
# explicit euler
return transientTerm.scalar(idx) + self.dt * self._discretizeSpatial(rhs)
return transientTerm.scalar + self.dt * self._discretizeSpatial(rhs)
else:
print(transientTerms)
raise NotImplementedError("Cannot discretize expression with more than one transient term")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment