[dune-pdelab] Initial conditions for instationary problems

["Jö Fahlke"]

Am Thu,  4. Jun 2015, 11:45:29 +0200 schrieb conf86 at web.de:
> Date: Thu, 4 Jun 2015 11:45:29 +0200
> From: conf86 at web.de
> To: Jö Fahlke <jorrit at jorrit.de>
> Cc: dune-pdelab at dune-project.org
> Subject: Re: [dune-pdelab] Initial conditions for instationary problems
> 
> > > I have another question regarding the dune-course examples.
> > > I'm trying to figure out where exactly the initial conditions for t=t_0 are set in example03. I see that in example03_Q1.hh there is a function calling the initial condition, but I don't really see the source.
> > > Or is it not set at all and 0 by default for this example? If so, is there an example-code how to set it?
> > 
> > I only have example03_Q2.hh or example02_Q1.hh available (because I'm still in
> > the train), but I guess the answer is as follows:
> 
> I'm sorry, I actually meant example03_Q2.hh.
> 
> > The dirichlet boundary conditions actually have to be consistent with the
> > initial conditions at t=0.  So in the examples we usually use the same
> > function for both: the stuff on *_bcextension.hh.  (BCExtension here means
> > "the function describing the boundary condition, extended to the interior of
> > the domain".)
> > 
> > The vector ok unknowns is then usually set to values interpolated from the
> > initial conditions by code like this:
> > 
> >   typedef BCExtension<GV,Real> G;                               // boundary value + extension
> >   G g(gv);
> >   Dune::PDELab::interpolate(g,gfs,u);                           // interpolate coefficient vector
> > 
> > The second line instantiates a BCExtension object "g", which in the third line
> > is interpoluted into the vector of unknowns "u", using the (grid) function
> > space "gfs".
> 
> Thank you for your answer 
> 
> I think, I understand the second part.
> 
> But how exactly can I set the initial values in bcextension.hh?
> My first guess would be something like this: 
> 
>  inline void evaluate(...){ 
>  ....
>  if (x[0]==0 && x[1]==0.5) y=2.0;
>   }
> 
> for an initial value u0(x,y)=2 for x=0, y=0.5. Is this correct?

Well, strictly speaking in a weird sense yes, but...

> So I guess I don't have to specify between initial values and boundary values and I could also write (still for a circle)? 
>  if (x[0]==0 && x[1]==1) y=1.0;

Your initial condition value must be defined on the whole domain[1].  You
probably have some analytic expression for it, in terms of the (global)
position x.  Just plug that in there.  I.e. if you wanted your initial
condition to be a Gaussian you would have something like

  y = std::exp(-x.two_norm()/(2*sigma*sigma))/sigma/std::sqrt(2*M_PI);

You don't need any "if" there to not set y in the outside -- evaluate will
simply never be called outside the domain.

As for the initial Dirichlet condition values: they _must_ match the initial
condition values at the boundary.  At least at t=0.  So as long as the
dirichlet condition values are independent of time, you can just reuse the
same function as for the initial conditions[2].

The way Dirichlet boundary conditions work is as follows.  The function in
*_bctypes.hh determines which DoFs are on the Dirichlet boundary and
constrains those DoFs.  That means that their values will not change when you
apply a linear solver.  And when you do a time step the value for the new
solution before the solve is copied from the old solution, and the solve will
not modify it so it will have the correct Dirichlet values.

Now it should be clear that the Dirichlet condition value is simply carried
along as you do time stepping.  But it needs to come from somewhere initially.
And that happens when you interpolate the initial condition, since at the
boundary the initial value and the dirichlet condition value must match.

I hope that helps,

Regards,
Jö.

[1] In theory.  In practice it is enough to have it defined at the
    interpolation points.  But you can't really predict those interpolation
    points.  So also in practice it is easiest to define the value on the
    whole domain.

[2] If you have time-dependent dirichlet boundary condition values, please get
    back to us, there are some more steps involved in that case.

-- 
Jorrit (Jö) Fahlke, Institute for Computational und Applied Mathematics,
University of Münster, Orleans-Ring 10, D-48149 Münster
Tel: +49 251 83 35146 Fax: +49 251 83 32729

In the beginning the Universe was created.  This has made a lot of
people very angry and been widely regarded as a bad move.
-- Douglas Adams
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 811 bytes
Desc: Digital signature
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20150608/e735a171/attachment.sig>