From c4f23fbe204ec54b440dfbd64e55d051e64584fd Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 9 Feb 2018 17:07:18 +0100
Subject: [PATCH] PressureTensor approach for phase field
---
phasefield/analytical.py | 126 ++++++++++++++++++++++-----------------
1 file changed, 70 insertions(+), 56 deletions(-)
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index c16028b6..d450b752 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -8,62 +8,6 @@ surfaceTensionSymbolName = "tau"
interfaceWidthSymbol = sp.Symbol("alpha")
-def functionalDerivative(functional, v, constants=None):
- """
- Computes functional derivative of functional with respect to v using Euler-Lagrange equation
-
- .. math ::
-
- \frac{\delta F}{\delta v} =
- \frac{\partial F}{\partial v} - \nabla \cdot \frac{\partial F}{\partial \nabla v}
-
- - assumes that gradients are represented by Diff() node (from Chapman Enskog module)
- - Diff(Diff(r)) represents the divergence of r
- - the constants parameter is a list with symbols not affected by the derivative. This is used for simplification
- of the derivative terms.
- """
- functional = expandUsingLinearity(functional, constants=constants)
- diffs = functional.atoms(Diff)
-
- diffV = Diff(v)
- assert diffV in diffs # not necessary in general, but for this use case this should be true
-
- nonDiffPart = functional.subs({d: 0 for d in diffs})
-
- partialF_partialV = sp.diff(nonDiffPart, v)
-
- dummy = sp.Dummy()
- partialF_partialGradV = functional.subs(diffV, dummy).diff(dummy).subs(dummy, diffV)
-
- result = partialF_partialV - Diff(partialF_partialGradV)
- return expandUsingLinearity(result, constants=constants)
-
-
-def coshIntegral(f, var):
- """Integrates a function f that has exactly one cosh term, from -oo to oo, by
- substituting a new helper variable for the cosh argument"""
- coshTerm = list(f.atoms(sp.cosh))
- assert len(coshTerm) == 1
- integral = sp.Integral(f, var)
- transformedInt = integral.transform(coshTerm[0].args[0], sp.Symbol("u", real=True))
- return sp.integrate(transformedInt.args[0], (transformedInt.args[1][0], -sp.oo, sp.oo))
-
-
-def discreteLaplace(field, index, dx):
- """Returns second order Laplace stencil"""
- dim = field.spatialDimensions
- count = 0
- result = 0
- for d in range(dim):
- for offset in (-1, 1):
- count += 1
- result += field.neighbor(d, offset)(index)
-
- result -= count * field.center(index)
- result /= dx ** 2
- return result
-
-
def symmetricSymbolicSurfaceTension(i, j):
"""Returns symbolic surface tension. The function is symmetric, i.e. interchanging i and j yields the same result.
If both phase indices i and j are chosen equal, zero is returned"""
@@ -181,6 +125,20 @@ def freeEnergyFunctionalNPhases(numPhases=None, surfaceTensions=symmetricSymboli
return result
+def separateIntoBulkAndInterface(freeEnergy):
+ """Separates the bulk and interface parts of a free energy
+
+ >>> F = freeEnergyFunctionalNPhases(3)
+ >>> bulk, inter = separateIntoBulkAndInterface(F)
+ >>> assert sp.expand(bulk - freeEnergyFunctionalNPhases(3, includeInterface=False)) == 0
+ >>> assert sp.expand(inter - freeEnergyFunctionalNPhases(3, includeBulk=False)) == 0
+ """
+ freeEnergy = freeEnergy.expand()
+ bulkPart = freeEnergy.subs({a: 0 for a in freeEnergy.atoms(Diff)})
+ interfacePart = freeEnergy - bulkPart
+ return bulkPart, interfacePart
+
+
def analyticInterfaceProfile(x, interfaceWidth=interfaceWidthSymbol):
"""Analytic expression for a 1D interface normal to x with given interface width
@@ -238,6 +196,62 @@ def createForceUpdateEquations(forceField, phiField, muField, dx=1):
return forceSweepEqs
+def functionalDerivative(functional, v, constants=None):
+ """
+ Computes functional derivative of functional with respect to v using Euler-Lagrange equation
+
+ .. math ::
+
+ \frac{\delta F}{\delta v} =
+ \frac{\partial F}{\partial v} - \nabla \cdot \frac{\partial F}{\partial \nabla v}
+
+ - assumes that gradients are represented by Diff() node (from Chapman Enskog module)
+ - Diff(Diff(r)) represents the divergence of r
+ - the constants parameter is a list with symbols not affected by the derivative. This is used for simplification
+ of the derivative terms.
+ """
+ functional = expandUsingLinearity(functional, constants=constants)
+ diffs = functional.atoms(Diff)
+
+ diffV = Diff(v)
+ #assert diffV in diffs # not necessary in general, but for this use case this should be true
+
+ nonDiffPart = functional.subs({d: 0 for d in diffs})
+
+ partialF_partialV = sp.diff(nonDiffPart, v)
+
+ dummy = sp.Dummy()
+ partialF_partialGradV = functional.subs(diffV, dummy).diff(dummy).subs(dummy, diffV)
+
+ result = partialF_partialV - Diff(partialF_partialGradV)
+ return expandUsingLinearity(result, constants=constants)
+
+
+def coshIntegral(f, var):
+ """Integrates a function f that has exactly one cosh term, from -oo to oo, by
+ substituting a new helper variable for the cosh argument"""
+ coshTerm = list(f.atoms(sp.cosh))
+ assert len(coshTerm) == 1
+ integral = sp.Integral(f, var)
+ transformedInt = integral.transform(coshTerm[0].args[0], sp.Symbol("u", real=True))
+ return sp.integrate(transformedInt.args[0], (transformedInt.args[1][0], -sp.oo, sp.oo))
+
+
+def discreteLaplace(field, index, dx):
+ """Returns second order Laplace stencil"""
+ dim = field.spatialDimensions
+ count = 0
+ result = 0
+ for d in range(dim):
+ for offset in (-1, 1):
+ count += 1
+ result += field.neighbor(d, offset)(index)
+
+ result -= count * field.center(index)
+ result /= dx ** 2
+ return result
+
+
def cahnHilliardFdEq(phaseIdx, phi, mu, velocity, mobility, dx, dt):
from pystencils.finitedifferences import transient, advection, diffusion, Discretization2ndOrder
cahnHilliard = transient(phi, phaseIdx) + advection(phi, velocity, phaseIdx) - diffusion(mu, mobility, phaseIdx)
--
GitLab