# Quasi-nonllinear Time-dependent Problem

Asked by Ted Kord on 2010-01-07

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:
2010-01-08
Last query:
2010-01-08
Last reply:
2010-01-08
 Nuno Lopes (ndlopes) said on 2010-01-07: #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.

 Nuno Lopes (ndlopes) said on 2010-01-07: #2

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

 Ted Kord (teddy-kord) said on 2010-01-07: #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.

 Nuno Lopes (ndlopes) said on 2010-01-07: #4

I can compile the form without any error...

 Ted Kord (teddy-kord) said on 2010-01-07: #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? Johan Hake (johan-hake) said on 2010-01-08: #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.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
>

 Ted Kord (teddy-kord) said on 2010-01-08: #7

Thanks Johan Hake, that solved my question.

To post a message you must log in.