dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Friday August 5 2011 13:26:52 Garth N. Wells wrote:
> On Aug 5 2011, Johan Hake wrote:
> >You are right. I will look at it this weekend.
> >
> > However, it "fixed" my problem with Epetra backend and Nonlinear problem
> > where I used the assignment operator when generating the Jacobian.
> 
> Can you post some simple code to test against?

Attached is a slightly modified version of the Cahn-hillard demo.

Johan
 
> Garth
> 
> >Johan
> >
> >On Aug 5, 2011, at 7:24, "Garth N. Wells" <gnw20@xxxxxxxxx> wrote:
> >> This doesn't look right. Looks like a shallow copy.
> >> 
> >> Garth
> >> 
> >> On 04/08/11 00:02, noreply@xxxxxxxxxxxxx wrote:
> >>> ------------------------------------------------------------
> >>> revno: 6069
> >>> committer: Johan Hake <hake.dev@xxxxxxxxx>
> >>> branch nick: work
> >>> timestamp: Wed 2011-08-03 20:58:24 -0700
> >>> 
> >>> message:
> >>>  Fix a bug in the assignment operator of EpetraMatrix
> >>> 
> >>> modified:
> >>>  dolfin/la/EpetraMatrix.cpp
> >>> 
> >>> --
> >>> lp:dolfin
> >>> https://code.launchpad.net/~dolfin-core/dolfin/main
> >>> 
> >>> Your team DOLFIN Core Team is subscribed to branch lp:dolfin. To
> >>> unsubscribe from this branch go to
> >>> https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription
> >> 
> >> _______________________________________________
> >> Mailing list: https://launchpad.net/~dolfin
> >> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> >> Unsubscribe : https://launchpad.net/~dolfin
> >> More help   : https://help.launchpad.net/ListHelp

"""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: 2009-06-20

# Begin demo

import random
from dolfin import *

parameters["linear_algebra_backend"] = "Epetra"

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self):
        random.seed(2 + MPI.process_number())
    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
        self.reset_sparsity = True
        self.A = Matrix()
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=self.A, reset_sparsity=self.reset_sparsity)
        A.assign(self.A)
        self.reset_sparsity = False

# 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 = UnitSquare(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("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 = 80*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()