DOLFIN

Quasi-nonllinear Time-dependent Problem

Asked by Ted Kord

Hi

I'm attempting to solve the following problem (using the C++ interface):

du/dt - div grad u = f(u)

Due to the non-linearity associated with f(u), I am using the parameters at time t to calculate it thus avoiding the use of a nonlinear root-finding algorithm.

The problem I have is how to define f(u).

Say, f(u) is:
      0.1/(1 + exp(u - 2) )

How do obtain the scalar value 'u' to pass as the argument to f(u) at every time step.

Thanks.

Ted

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Nuno Lopes (ndlopes) said : 		#1

I guess it could be done in the ufl form file:

Something like (appart from the time discretization algorithm)

element = FiniteElement("Lagrange", triangle, 1)

v = TestFunction(element)
u = TrialFunction(element) # Unkown Function

u_old =Function(element) ##Known Function from the last step
##

a=u*v*dx+inner(grad(u),grad(v))*dx
L=u_old*v*dx+v*0.1/(exp(u_old-2))*dx

Hope it helps.

Revision history for this message
Nuno Lopes (ndlopes) said : 		#2

I forgot to multiply the timestep in the forms...

Revision history for this message
Ted Kord (teddy-kord) said : 		#3

That was my initial approach but I get the error:

L=u_old*v*dx+v*0.1/(exp(u_old-2))*dx
"TypeError: a float is required"

That's why I was wondering how to get to the scalar value (in the ufl form or in a subclass of Expression) assuming it can't be done a different way.

Revision history for this message
Nuno Lopes (ndlopes) said : 		#4

I can compile the form without any error...

Revision history for this message
Ted Kord (teddy-kord) said : 		#5

It seems the form compiler wasn't too happy with the exp(...) so I replace it with e**(...) .

I no longer get the Typeerror but now I get this at Compiler Stage 4:

Compiler stage 4: Generating code
---------------------------------
  Generating code using quadrature representation.

  Generating code for cell integrals.
  Transforming UFL integrand...power does not support this exponent: Division(Sum(IntValue(25, (), (), {}), Product(IntValue(-1, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), FloatValue(5.9800000000000004, (), (), {}))

*** FFC: power does not support this exponent: Division(Sum(IntValue(25, (), (), {}), Product(IntValue(-1, (), (), {}), Coefficient(FiniteElement('Lagrange', Cell('triangle', 1, Space(2)), 1), 0))), FloatValue(5.9800000000000004, (), (), {}))

Is there a workaround for this?

Revision history for this message
Best Johan Hake (johan-hake) said : 		#6

On Thursday 07 January 2010 09:18:00 Ted Kord wrote:
> New question #96571 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/96571
>
> Hi
>
> I'm attempting to solve the following problem (using the C++ interface):
>
> du/dt - div grad u = f(u)
>
> Due to the non-linearity associated with f(u), I am using the parameters at
> time t to calculate it thus avoiding the use of a nonlinear root-finding
> algorithm.
>
> The problem I have is how to define f(u).
>
> Say, f(u) is:
> 0.1/(1 + exp(u - 2) )
>
> How do obtain the scalar value 'u' to pass as the argument to f(u) at every
> time step.

Have you tried to just define a scalar Function or Coefficient (new name) in
ufl for f(). Then you subclass Expression in C++, something like:

  class F: public Expression
  {
    Function * _u;
  public:
    F:Expression(Function& u):Expression(),_u(&u){}

    void eval(Array<double>& values, const Array<const double>& x) const
    {
      double val;
      _u->eval(val, x);
      values[0] = 0.1/(1+exp(val-2));
    }
  };

  ...

  Function u_0(V);
  F f(u_0);

You should then have the values of the previous function available.

Johan

> Thanks.
>
> Ted
>

Revision history for this message
Ted Kord (teddy-kord) said : 		#7

Thanks Johan Hake, that solved my question.