dolfin team mailing list archive

[jesperc <jesperc@xxxxxx>]

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;

}