Skip to content
Snippets Groups Projects
Commit 1e5fc1b9 authored by Martin Bauer's avatar Martin Bauer
Browse files

Renamed misleading discretization function

parent e475d0e1
No related branches found
No related tags found
No related merge requests found
...@@ -61,7 +61,8 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt ...@@ -61,7 +61,8 @@ def freeEnergyFunction3Phases(orderParameters=None, interfaceWidth=interfaceWidt
def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymbolicSurfaceTension, def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymbolicSurfaceTension,
interfaceWidth=interfaceWidthSymbol, orderParameters=None, interfaceWidth=interfaceWidthSymbol, orderParameters=None,
includeBulk=True, includeInterface=True, symbolicLambda=False): includeBulk=True, includeInterface=True, symbolicLambda=False,
symbolicDependentVariable=False):
r""" r"""
Returns a symbolic expression for the free energy of a system with N phases and Returns a symbolic expression for the free energy of a system with N phases and
specified surface tensions. The total free energy is the sum of a bulk and an interface component. specified surface tensions. The total free energy is the sum of a bulk and an interface component.
...@@ -89,6 +90,9 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli ...@@ -89,6 +90,9 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
Parameter useful for viewing / debugging the function Parameter useful for viewing / debugging the function
:param includeBulk: if false no bulk term is added :param includeBulk: if false no bulk term is added
:param includeInterface:if false no interface contribution is added :param includeInterface:if false no interface contribution is added
:param symbolicLambda: surface energy coefficient is represented by symbol, not in expanded form
:param symbolicDependentVariable: last phase variable is defined as 1-otherPhaseVars, if this is set to True
it is represented by phi_A for better readability
""" """
assert not (numPhases is None and orderParameters is None) assert not (numPhases is None and orderParameters is None)
if orderParameters is None: if orderParameters is None:
...@@ -97,7 +101,10 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli ...@@ -97,7 +101,10 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
phi = orderParameters phi = orderParameters
numPhases = len(phi) + 1 numPhases = len(phi) + 1
phi = tuple(phi) + (1 - sum(phi),) if not symbolicDependentVariable:
phi = tuple(phi) + (1 - sum(phi),)
else:
phi = tuple(phi) + (sp.Symbol("phi_D"), )
# Compared to handwritten notes we scale the interface width parameter here to obtain the correct # 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 # equations for the interface profile and the surface tensions i.e. to pass tests
...@@ -109,11 +116,11 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli ...@@ -109,11 +116,11 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
def lambdaCoeff(k, l): def lambdaCoeff(k, l):
if symbolicLambda: if symbolicLambda:
return sp.Symbol("Lambda_%d%d" % ((k,l) if k < l else (l,k))) return sp.Symbol("Lambda_%d%d" % ((k, l) if k < l else (l, k)))
N = numPhases - 1 N = numPhases - 1
if k == l: if k == l:
assert surfaceTensions(l, l) == 0 assert surfaceTensions(l, l) == 0
return 3 / sp.sqrt(2) * interfaceWidth * (surfaceTensions(k, N) + surfaceTensions(l, N) - surfaceTensions(k, l)) return 3 / sp.sqrt(2) * interfaceWidth * (surfaceTensions(k, N) + surfaceTensions(l, N) - surfaceTensions(k, l))
def bulkTerm(i, j): def bulkTerm(i, j):
return surfaceTensions(i, j) / 2 * (f(phi[i]) + f(phi[j]) - f(phi[i] + phi[j])) return surfaceTensions(i, j) / 2 * (f(phi[i]) + f(phi[j]) - f(phi[i] + phi[j]))
...@@ -234,7 +241,7 @@ def coshIntegral(f, var): ...@@ -234,7 +241,7 @@ def coshIntegral(f, var):
return sp.integrate(transformedInt.args[0], (transformedInt.args[1][0], -sp.oo, sp.oo)) return sp.integrate(transformedInt.args[0], (transformedInt.args[1][0], -sp.oo, sp.oo))
def discretizeSecondDerivatives(term, dx=1): def finiteDifferences2ndOrder(term, dx=1):
"""Substitutes symbolic integral of field access by second order accurate finite differences. """Substitutes symbolic integral of field access by second order accurate finite differences.
The only valid argument of Diff objects are field accesses (usually center field accesses)""" The only valid argument of Diff objects are field accesses (usually center field accesses)"""
def diffOrder(e): def diffOrder(e):
......
import sympy as sp import sympy as sp
from lbmpy.phasefield.analytical import chemicalPotentialsFromFreeEnergy, substituteLaplacianBySum, \ from lbmpy.phasefield.analytical import chemicalPotentialsFromFreeEnergy, substituteLaplacianBySum, \
discretizeSecondDerivatives, forceFromPhiAndMu, symmetricTensorLinearization, pressureTensorFromFreeEnergy, \ finiteDifferences2ndOrder, forceFromPhiAndMu, symmetricTensorLinearization, pressureTensorFromFreeEnergy, \
forceFromPressureTensor forceFromPressureTensor
...@@ -14,7 +14,7 @@ def muKernel(freeEnergy, orderParameters, phiField, muField, dx=1): ...@@ -14,7 +14,7 @@ def muKernel(freeEnergy, orderParameters, phiField, muField, dx=1):
chemicalPotential = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters) chemicalPotential = chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters)
chemicalPotential = substituteLaplacianBySum(chemicalPotential, dim) chemicalPotential = substituteLaplacianBySum(chemicalPotential, dim)
chemicalPotential = chemicalPotential.subs({op: phiField(i) for i, op in enumerate(orderParameters)}) chemicalPotential = chemicalPotential.subs({op: phiField(i) for i, op in enumerate(orderParameters)})
return [sp.Eq(muField(i), discretizeSecondDerivatives(mu_i, dx)) for i, mu_i in enumerate(chemicalPotential)] return [sp.Eq(muField(i), finiteDifferences2ndOrder(mu_i, dx)) for i, mu_i in enumerate(chemicalPotential)]
def forceKernelUsingMu(forceField, phiField, muField, dx=1): def forceKernelUsingMu(forceField, phiField, muField, dx=1):
...@@ -22,7 +22,7 @@ def forceKernelUsingMu(forceField, phiField, muField, dx=1): ...@@ -22,7 +22,7 @@ def forceKernelUsingMu(forceField, phiField, muField, dx=1):
assert muField.indexDimensions == 1 assert muField.indexDimensions == 1
force = forceFromPhiAndMu(phiField.vecCenter, mu=muField.vecCenter, dim=muField.spatialDimensions) force = forceFromPhiAndMu(phiField.vecCenter, mu=muField.vecCenter, dim=muField.spatialDimensions)
return [sp.Eq(forceField(i), return [sp.Eq(forceField(i),
discretizeSecondDerivatives(f_i, dx)).expand() for i, f_i in enumerate(force)] finiteDifferences2ndOrder(f_i, dx)).expand() for i, f_i in enumerate(force)]
def pressureTensorKernel(freeEnergy, orderParameters, phiField, pressureTensorField, dx=1): def pressureTensorKernel(freeEnergy, orderParameters, phiField, pressureTensorField, dx=1):
...@@ -34,7 +34,7 @@ def pressureTensorKernel(freeEnergy, orderParameters, phiField, pressureTensorFi ...@@ -34,7 +34,7 @@ def pressureTensorKernel(freeEnergy, orderParameters, phiField, pressureTensorFi
eqs = [] eqs = []
for index, linIndex in indexMap.items(): for index, linIndex in indexMap.items():
eq = sp.Eq(pressureTensorField(linIndex), eq = sp.Eq(pressureTensorField(linIndex),
discretizeSecondDerivatives(p[index], dx).expand()) finiteDifferences2ndOrder(p[index], dx).expand())
eqs.append(eq) eqs.append(eq)
return eqs return eqs
...@@ -45,7 +45,7 @@ def forceKernelUsingPressureTensor(forceField, pressureTensorField, dx=1): ...@@ -45,7 +45,7 @@ def forceKernelUsingPressureTensor(forceField, pressureTensorField, dx=1):
p = sp.Matrix(dim, dim, lambda i, j: pressureTensorField(indexMap[i,j] if i < j else indexMap[j, i])) p = sp.Matrix(dim, dim, lambda i, j: pressureTensorField(indexMap[i,j] if i < j else indexMap[j, i]))
f = forceFromPressureTensor(p) f = forceFromPressureTensor(p)
return [sp.Eq(forceField(i), discretizeSecondDerivatives(f_i, dx).expand()) return [sp.Eq(forceField(i), finiteDifferences2ndOrder(f_i, dx).expand())
for i, f_i in enumerate(f)] for i, f_i in enumerate(f)]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment