Here are the two codes for FEniCS (demo_cahn-hilliard.py) and FiPy
(spinodal2D.py). I have commented out the file saving part in both the programs
so that we can only compare the computation time. The mesh sizes and time
stepping are also equivalent. I still believe that there is something I am
doing wrong, so please have a careful look. Thanks a lot.
"""This demo illustrates how to use of DOLFIN for solving the Cahn-Hilliard
equation, which is a time-dependent nonlinear PDE """
# Copyright (C) 2009 Garth N. Wells
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# First added: 2009-06-20
# Last changed: 2013-11-20
# Begin demo
import time
import random
from dolfin import *
# Class representing the intial conditions
class InitialConditions(Expression):
def __init__(self):
random.seed(2 + MPI.rank(mpi_comm_world()))
def eval(self, values, x):
values[0] = 0.63 + 0.02*(0.5 - random.random())
values[1] = 0.0
def value_shape(self):
return (2,)
# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
# Model parameters
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"
# Create mesh and define function spaces
mesh = UnitSquareMesh(96, 96)
V = FunctionSpace(mesh, "Lagrange", 1)
ME = V*V
# Define trial and test functions
du = TrialFunction(ME)
q, v = TestFunctions(ME)
# Define functions
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
# Create intial conditions and interpolate
u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)
# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
#file = File("output.pvd", "compressed")
# Step in time
t = 0.0
T = 50*dt
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
#file << (u.split()[0], t)
#plot(u.split()[0])
#interactive()
#!/usr/bin/env python
from fipy.tools.parser import parse
numberOfElements = parse('--numberOfElements', action = 'store', type = 'int', default = 9216) # 96*96
numberOfSteps = parse('--numberOfSteps', action = 'store', type = 'int', default = 50)
import fipy.tools.numerix as numerix
nx = int(numerix.sqrt(numberOfElements))
ny = int(numerix.sqrt(numberOfElements))
steps = numberOfSteps
dx = 1.0 ##/nx
dy = 1.0 ##/ny
asq = 1.0
epsilon = 1.0 ##0.5e-2
diffusionCoeff = 1 ## e-5
from fipy.meshes.grid2D import Grid2D
mesh = Grid2D(dx = dx, dy = dy, nx = nx, ny = ny)
from fipy.variables.cellVariable import CellVariable
from fipy.tools.numerix import random
var = CellVariable(name = "phase field",
mesh = mesh,
value = random.random(nx * ny))
faceVar = var.getArithmeticFaceValue()
doubleWellDerivative = asq * ( 1 - 6 * faceVar * (1 - faceVar))
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.transientTerm import TransientTerm
diffTerm2 = ImplicitDiffusionTerm(coeff = (diffusionCoeff * doubleWellDerivative,))
diffTerm4 = ImplicitDiffusionTerm(coeff = (diffusionCoeff, -epsilon**2))
eqch = TransientTerm() - diffTerm2 - diffTerm4
from fipy.solvers.pysparse.linearPCGSolver import LinearPCGSolver
from fipy.solvers.pysparse.linearLUSolver import LinearLUSolver
##solver = LinearLUSolver(tolerance = 1e-15,steps = 1000)
solver = LinearLUSolver(tolerance = 1e-10,iterations = 1000)
from fipy.boundaryConditions.fixedValue import FixedValue
from fipy.boundaryConditions.fixedFlux import FixedFlux
from fipy.boundaryConditions.nthOrderBoundaryCondition import NthOrderBoundaryCondition
BCs = (FixedFlux(mesh.getFacesRight(), 0),
FixedFlux(mesh.getFacesLeft(), 0),
NthOrderBoundaryCondition(mesh.getFacesLeft(), 0, 3),
NthOrderBoundaryCondition(mesh.getFacesRight(), 0, 3),
NthOrderBoundaryCondition(mesh.getFacesTop(), 0, 3),
NthOrderBoundaryCondition(mesh.getFacesBottom(), 0, 3))
if __name__ == '__main__':
import time
start = time.time()
#import fipy.viewers
#viewer = fipy.viewers.make(vars = var, limits = {'datamin': 0., 'datamax': 1.0})
#viewer.plot()
#dexp=-5
dt = 5.0e-6
for step in range(steps):
#dt = numerix.exp(dexp)
#dt = min(100, dt)
#dexp += 0.01
#var.updateOld()
eqch.sweep(var, boundaryConditions = BCs, solver = solver, dt = dt)
end = time.time()
print "Time elapsed:", end - start
'''
if __name__ == '__main__':
viewer.plot()
print 'step',step,'dt',dt
'''
def _run():
pass
On Aug 13, 2014, at 2:32 AM, Jan Blechta <[email protected]> wrote:
> Could you provide the codes? It would be interesting...
>
> Jan
>
>
> On Tue, 12 Aug 2014 19:59:45 -0400
> Aniruddha Jana <[email protected]> wrote:
>
>> Thanks for the detailed thoughts, that helps. For the time stepping,
>> I used the same time stepping for FiPy as given in the Fenics C-H
>> demo (a constant increment). I made the two codes as similar as
>> possible from the front end before the comparison.
>>
>>
>> On Aug 12, 2014, at 3:30 PM, Mike Welland <[email protected]>
>> wrote:
>>
>>> One contributing factor could be the difference between finite
>>> element vs. finite volume. In the fenics CH demo, the 4th order
>>> diff eq is split into 2nd order eqs. for reasons discussed under
>>> section 5.1.1. Mixed Form on the doc page. FiPy's demo doesn't
>>> seem to do that (at least as far as I can see). Finite difference
>>> codes also don't need to split.
>>>
>>> Now if that could explain a factor of almost 40 is another matter
>>> entirely. You would have to dig into things like time stepping,
>>> error control, etc. .e.g: It looks like FiPy uses an exponential
>>> time step whereas the fenics version uses a constant.
>>>
>>> Depending on what you want to do ultimately, bear in mind issues
>>> like parallelization, mesh refinement, supported linear backend
>>> (fenics = PETSc, dunno about FiPy) etc.
>>>
>>> Mike
>>>
>>>
>>> On Tue, Aug 12, 2014 at 2:13 PM, Aniruddha Jana <[email protected]>
>>> wrote: I ran with same mesh sizes and for equal time steps.
>>>
>>> On August 12, 2014 1:57:30 PM EDT, Jan Blechta
>>> <[email protected]> wrote: If FyPi Canh-Hilliard example
>>> is this one
>>> http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html
>>>
>>>
>>> then the reason is very simple:
>>>
>>> # FEniCS
>>> mesh = UnitSquareMesh(96, 96)
>>>
>>> # FiPy
>>> __name__ == "__main__":
>>> nx = ny = 20
>>> else:
>>> nx = ny = 10
>>> mesh = Grid2D(nx=nx, ny=ny, dx=0.25, dy=0.25)
>>>
>>>
>>>
>>> Also the function space may be different if FiPy's 'CellVariable' is
>>> something-like piece-wise constants.
>>>
>>> Jan
>>>
>>>
>>> On Tue, 12 Aug 2014 10:42:50 -0400
>>> Aniruddha Jana <[email protected]> wrote:
>>>
>>>
>>>
>>> Hello,
>>>
>>> I am trying to learn FEniCS, and have been using FiPy so far. I !
>>> ran
>>>
>>> the Python Cahn-Hilliard example. The program took around 80
>>> seconds to run in serial, while a FiPy Cahn-Hillard program with
>>> similar size and settings took only 2.72 seconds. I think I am
>>> making some mistake
>>>
>>>
>>> here as I expected FEniCS to be better than FiPy in terms of
>>> speed.
>>> Can somebody please comment on the speed and memory issues,
>>> especially in comparison to FiPy? Since I am trying to learn using
>>> FEniCS, I would appreciate any such comments.
>>>
>>>
>>>
>>> Many thanks,
>>> Aniruddha
>>>
>>> fenics mailing list
>>> [email protected]
>>> http://fenicsproject.org/mailman/listinfo/fenics
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> fenics mailing list
>>> [email protected]
>>> http://fenicsproject.org/mailman/listinfo/fenics
>>>
>>>
>>
>
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
