Unstable solution

[Krzysztof Bzowski]

Hi,

I have a problem with obtaining stable results of the time dependent problem. In the beginning results are correct but after few iteration solution start changing unexpectedly. Maybe somebody had similar problem and found solution.

I am trying to solve Fick's Second law:

Using Galerkin's method I have obtain following formulas:


T0 is solution vector obtained from previous time step. T1 is current unknow.

My assemble function is based on a solution of Laplace problem (from deal.ii examples) and looks like: http://pastebin.com/Ls5nS95p

Mesh was obtained from gmsh and it is pretty dense. I have used second order elements. I have applied Dirichlet boundary condition for several elements in the corners using VectorTools::interpolate_boundary_values, the value of the boundary condition is changing with each time step.

Values calculated in previous time steps are copied and used as a T0 (initial condition) in current time step.

At the beginning results look stable. But after few iterations solution starts changing unpredictable: http://imgur.com/a/oDvmP

Code is exactly the same. Each iteration step looks identical.

Using denser grid cause that solution looks good till the end of simulation. But denser grid is connected with a longer computation time which I would like to avoid. It's strange for me that errors occurs after some time. I've got only basic FEM training so maybe I miss something. Can I fight somehow with that errors?

Best regards,
Krzysztof

[Wolfgang Bangerth]

Krzysztof,
I apparently can't inline your formulas, but I do not understand how your time
stepping scheme works. You list this:

(2H + 3/dt) T1 = -(H + 3/dt C) T0

First, there is surely a C missing after 3/dt in the parentheses on the left.
Secondly, I don't understand where this formula is coming from. I can reorder
it to

C (T1 + T0)/dt + 1/3 H (2 T1 + T0) = 0

What kind of time discretization is this supposed to be? For sure, the first
term does not represent a time derivative (it should be T1-T0 instead of T1+T0).


Other than that, there is a description of how I debug time dependent problems
here:

http://www.dealii.org/developer/doxygen/deal.II/step_26.html#Verifyingwhetherthecodeiscorrect

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Krzysztof Bzowski]

On date 2013-07-17 15:57, Wolfgang Bangerth wrote:
>
> Krzysztof,
> I apparently can't inline your formulas, but I do not understand how
> your time stepping scheme works. You list this:
>
> (2H + 3/dt) T1 = -(H + 3/dt C) T0
>
> First, there is surely a C missing after 3/dt in the parentheses on the
> left. Secondly, I don't understand where this formula is coming from. I
> can reorder it to
>
> C (T1 + T0)/dt + 1/3 H (2 T1 + T0) = 0
>
> What kind of time discretization is this supposed to be? For sure, the
> first term does not represent a time derivative (it should be T1-T0
> instead of T1+T0).
>

I copied that formulas from my notes and of course I made an error.
It should be:
(2H + (3/dt)*C)*T1 + (H - (3/dt)*C)*T0 + 3P = 0
And after reordered it gives:

C(T1-T0)/dt + H(T0/3 + 2/3*T1) = 0
P contains external fluxes so it is equal 0 in my case.

The formula is based on the method of weighted residuals.

I have checked - I have implemented correct formula in my code. Do you
think it is wrong?



Best regards,
Krzysztof Bzowski

[Wolfgang Bangerth]

> C(T1-T0)/dt + H(T0/3 + 2/3*T1) = 0
> P contains external fluxes so it is equal 0 in my case.
>
> The formula is based on the method of weighted residuals.
>
> I have checked - I have implemented correct formula in my code. Do you think
> it is wrong?

This may be correct. Can you try the Crank-Nicolson method, which is

C(T1-T0)/dt + H(T0/2 + T1/2) = 0

? If that still yields problems, then there simply must be a bug in your program.

[Krzysztof Bzowski]

On date 2013-07-17 17:25, Wolfgang Bangerth wrote:
> This may be correct. Can you try the Crank-Nicolson method, which is
>
> C(T1-T0)/dt + H(T0/2 + T1/2) = 0
>
> ? If that still yields problems, then there simply must be a bug in your
> program.
>
> W.
>

Thank you for you time and help. I have checked Crank-Nicolson method
but errors still occurred.

I was extensively debugging my application and found possible cause.
It was diffusion coefficient, which value depends on time.

In the case where the value of the coefficient is constant - everything
went smoothly without errors. But when I changed it to time dependent -
the problem occurs.
Physics behind my simulation requires a time dependent diffusion
coefficient (it needs decrease with time).

Best regards,
Krzysztof