dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

What is the problem with the second option? Do you get an error
message or does it simply not converge?

/Anders


On Mon, Jun 04, 2007 at 01:54:38PM +0200, jesperc wrote:
> Hi everyone,
> 
> I have a problem with the non-linear solver. I want to solve a 
> non-linear pde of the form f(u,grad(u))*dot(grad(u),grad(v)) = 0 for all 
> v. If I define f inside the form-file everything works fine, but if I 
> define f as a function which is calculated outside the form-file it only 
> works when f=f(u) but not when f=f(grad(u)).
> 
> I have attached an example of this discrepancy.
> 
> Regards, Jesper
> 
> 

> # Copyright (C) 2007 Jesper Carlsson.
> # Licensed under the GNU GPL Version 2.
> #
> # Modified by Jesper Carlsson
> # First added:  2007-04-03
> # Last changed: 2007-04-03
> 
> # The linear form L(v) for the gradient product at each cell.
> #
> # Compile this form with FFC: ffc MassMat.form.
> 
> element0 = FiniteElement("Discontinuous Lagrange", "triangle", 0)
> element1 = FiniteElement("Lagrange", "triangle", 1)
> 
> v = TestFunction(element0)
> f = Function(element1)
> 
> # Used to assemble gradients of functions
> L = dot(grad(f), grad(f))*v*dx          

> # Copyright (C) 2007 Jesper Carlsson.
> # Licensed under the GNU GPL Version 2.
> #
> # Modified by Jesper Carlsson
> # First added:  2007-05-25
> # Last changed: 2007-06-01
> 
> element = FiniteElement("Lagrange", "triangle", 1)
> 
> v = TestFunction(element)
> u = TrialFunction(element)
> 
> # Alternative 1 (works)
> 
> u0= Function(element)
> 
> a = (1+u0*u0)*dot(grad(u),grad(v))*dx + 2*u0*u*dot(grad(u0),grad(v))*dx
> L = -(1+u0*u0)*dot(grad(u0),grad(v))*dx

> # Copyright (C) 2007 Jesper Carlsson.
> # Licensed under the GNU GPL Version 2.
> #
> # Modified by Jesper Carlsson
> # First added:  2007-05-25
> # Last changed: 2007-06-01
> 
> element = FiniteElement("Lagrange", "triangle", 1)
> 
> v = TestFunction(element)
> u = TrialFunction(element)
> 
> # Alternative 2 (?)
> 
> u0= Function(element)
> f = Function(element)  # f=1+u0^2
> fp = Function(element) # fp=1
> 
> a = f*dot(grad(u),grad(v))*dx + 2*fp*u0*u*dot(grad(u0),grad(v))*dx
> L = -f*dot(grad(u0),grad(v))*dx

> # Copyright (C) 2007 Jesper Carlsson.
> # Licensed under the GNU GPL Version 2.
> #
> # Modified by Jesper Carlsson
> # First added:  2007-05-25
> # Last changed: 2007-06-01
> 
> element = FiniteElement("Lagrange", "triangle", 1)
> 
> v = TestFunction(element)
> u = TrialFunction(element)
> 
> # Alternative 3 (works)
> 
> u0= Function(element)
> 
> a = (1+dot(grad(u0),grad(u0)))*dot(grad(u),grad(v))*dx + 2*dot(grad(u0),grad(u))*dot(grad(u0),grad(v))*dx
> L = -(1+dot(grad(u0),grad(u0)))*dot(grad(u0),grad(v))*dx
> 
> 

> # Copyright (C) 2007 Jesper Carlsson.
> # Licensed under the GNU GPL Version 2.
> #
> # Modified by Jesper Carlsson
> # First added:  2007-05-25
> # Last changed: 2007-06-01
> 
> element0 = FiniteElement("Discontinuous Lagrange", "triangle", 0)
> element1 = FiniteElement("Lagrange", "triangle", 1)
> 
> v = TestFunction(element1)
> u = TrialFunction(element1)
> 
> # Alternative 4 (does not work)
> 
> u0= Function(element1)
> f= Function(element0)  # f = 1 + |grad(u)|^2
> fp= Function(element0) # fp = 1
> 
> a = f*dot(grad(u),grad(v))*dx + 2*fp*dot(grad(u0),grad(u))*dot(grad(u0),grad(v))*dx
> L = -f*dot(grad(u0),grad(v))*dx

> // Copyright (C) 2007 Jesper Carlsson.
> // Licensed under the GNU GPL Version 2.
> //
> // Modified by Jesper Carlsson
> // First added:  2007-05-29
> // Last changed: 2007-05-29
> //
> 
> #include <dolfin.h>
> #include "dolfin/Test1.h"
> #include "dolfin/Test2.h"
> #include "dolfin/Test3.h"
> #include "dolfin/Test4.h"
> #include "dolfin/GradProd.h"
>   
> using namespace dolfin;
> 
> int TEST = 4;
> 
> // Dirichlet boundary condition
> class DirichletBC : public BoundaryCondition {
>   void eval(BoundaryValue& value, const Point& p, unsigned int i) {
>     if ( fabs(p.x()) < DOLFIN_EPS )
>       value = 0.25-pow(p.y()-0.5, 2);
>     else
>       value = 0.0;
>   }
> };
> 
> // User defined nonlinear problem 
> class MyNonlinearProblem : public NonlinearProblem
> {
> public:
>   
>   // Constructor 
>   MyNonlinearProblem(Mesh& mesh, BoundaryCondition& dbc, Function& U,
>                      Function& F, Function& Fprim) 
>     : NonlinearProblem(), _mesh(&mesh), _dbc(&dbc), _U(&U), _F(&F), _Fprim(&Fprim)
>   {
> 
>     // Create forms
> 
>     if (TEST==1)
>       {
>         // Alternative 1 (OK)
>         a = new Test1::BilinearForm(U);
>         L = new Test1::LinearForm(U);
>       }
>     else if (TEST==2)
>       {
>         // Alternative 2 (OK)
>         a = new Test2::BilinearForm(U, F, Fprim);
>         L = new Test2::LinearForm(U, F);
>       }
>     else if (TEST==3)
>       {
>         // Alternative 1 (OK)
>         a = new Test3::BilinearForm(U);
>         L = new Test3::LinearForm(U);
>       }
>     else if (TEST==4)
>       {
>         // Alternative 4 (NOT OK)
>         a = new Test4::BilinearForm(U, F, Fprim);
>         L = new Test4::LinearForm(U, F);
>       }
>     
>     gp = new GradProd::LinearForm(U);
>     
>     // Initialise solution vector u
>     U.init(mesh, a->test());
> 
>     // Initialise cell-wise defined vectors
>     F.init(mesh, gp->test());
>     Fprim.init(mesh, gp->test());
> 
>   }
>     
>   // Destructor 
>   ~MyNonlinearProblem()
>   {
>     delete a;
>     delete L;
>     delete gp;
>   }
>   
>   // User defined assemble of Jacobian and residual vector 
>   void form(GenericMatrix& A, GenericVector& b, const GenericVector& x)
>   {
> 
>     // Assemble gradients
>     dolfin_log(false);
>     FEM::assemble(*gp, grad2, *_mesh);
>     dolfin_log(true);
> 
>     // Calculate f and fp
>     dolfin_assert(_F->vector().size()==grad2.size());
>     dolfin_assert(_Fprim->vector().size()==grad2.size());
> 
>     for (int i=0; i<grad2.size(); ++i)
>       {
>         if (TEST==2)
>           {
>             // Alternative 2
>             _F->vector()(i) = 1+pow(_U->vector().get(i), 2);
>             _Fprim->vector()(i) = 1;
>           }
>         else if (TEST==4)
>           {
>             // Alternative 4
>             _F->vector()(i) = 1+grad2(i);
>             _Fprim->vector()(i) = 1;
>           }
>         
>       }
> 
>     dolfin_log(false);
>     FEM::assemble(*a, *L, A, b, *_mesh);
>     FEM::applyBC(A, *_mesh, a->test(), *_dbc);
>     FEM::applyResidualBC(b, x, *_mesh, a->test(), *_dbc);
>     dolfin_log(true);
>   }
>   
> private:
>   
>   // Pointers to forms, mesh and boundary conditions
>   BilinearForm *a;
>   LinearForm *L, *gp;
>   Mesh* _mesh;
>   BoundaryCondition* _dbc;
>   Function *_U, *_F, *_Fprim;
>   Vector grad2;
>   
> };
> 
> int main()
> {
> 
>   // Set up problem
>   UnitSquare mesh(64, 64);
>   DirichletBC dbc;
>   Function U, F, Fprim;
>   
>   // Create user-defined nonlinear problem
>   MyNonlinearProblem nonlinear_problem(mesh, dbc, U, F, Fprim);
> 
>   // Create nonlinear solver and set parameters
>   NewtonSolver nonlinear_solver;
>   nonlinear_solver.set("Newton maximum iterations", 50);
>   nonlinear_solver.set("Newton relative tolerance", 1e-10);
>   nonlinear_solver.set("Newton absolute tolerance", 1e-10);
> 
>   
>   // Solve nonlinear problem
>   Vector& u = U.vector();
>   nonlinear_solver.solve(nonlinear_problem, u);
>     
>   // Save function to file
>   File file("nonlinear_poisson.pvd");
>   file << U;
>   
>   return 0;
> 
> }

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev