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" <gn...@cam.ac.uk> wrote: > >> This doesn't look right. Looks like a shallow copy. > >> > >> Garth > >> > >> On 04/08/11 00:02, nore...@launchpad.net wrote: > >>> ------------------------------------------------------------ > >>> revno: 6069 > >>> committer: Johan Hake <hake....@gmail.com> > >>> 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@lists.launchpad.net > >> 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()
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
