dolfin team mailing list archive

["N.A. Borggren" <question159192@xxxxxxxxxxxxxxxxxxxxx>]

Question #159192 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/159192

N.A. Borggren proposed the following answer:
Hello,
  I've been working on related Fokker Planck type problems (the 'adjoint' of
what you wrote down) and have encountered this same phenomena.  I have not
managed to find a way to get directly the steady state but perhaps you can
learn from what I've tried.  The fastest way I've been able to converge to a
steady state is by still solving in the time dependent way (like the
advection/diffusion demo) but using an initial condition 'as close as
possible' to the steady state.  For example obtaining an initial guess by
histogramming a few of the corresponding dynamical system or langevin
trajectories and solving poisson's equation gives an initial distribution
that converges quite rapidly.  Also artificially increasing the diffusion
(high temperature limit) can speed up the simulation and then you can turn
it back down to the value of interest when it seems closer to converged.
  I know those aren't the answers you were hoping for, but that is the best
I've managed to work out so far.  Part of the issue, I believe, is the
'steady state' is, not necessarily the same for all initial conditions.  It
is often the case that multiple initial conditions asymptotically
convergence to the same steady state but that technically takes infinitely
long.   The equation is still a linear equation, a unique initial condition
gives a unique final condition, and two different initial conditions that
seem 'converged' will actual differ a bit in the details so some tolerance
would need to be defined and set to consider these two solutions identical.
  Also if the boundary conditions are at all absorbing then the steady state
actually is the trivial solution unless there is an explicit driving force.
(take t to infinity in, say, Equation 2.40 Zwanzig, Nonequilibrium
statistical mechanics).  So depending on your boundary condition dolfin and
your form might actually be working properly.
  I was hoping to avoid all that with a reflecting boundary condition
case, (A_i
* u)n_i + (d/dx_j (B^2_ij * u)n_i=0, on the mesh boundary with n_i the
normal to the mesh.  This is how I keep the probability from going to zero
in the time dependent case as well since it conserves the probability inside
the mesh domain.  For the steady state I have also tried to at least avoid
the trivial solution by forcing a normalization condition of 1 with a
lagrange multiplier but haven't yet succeeded with that approach either.
The iterative method of finding a good guess from the dynamics for the
initial condition and propagating that in time has been the fastest, albeit
slightly unsatisfying method for me.
  I think for these steady states, as well as the time dependent case,
Fokker-Planck equations would work a lot better if one were able to take
into account that we need only be concerned with positive functions and
could restrict the finite elements used to insist on positivity of a
function.  Are there strictly positive finite elements available or good
ways to force the standard elements to stay positive?  It could be very
useful here.  This has also been the hold up I've had in applying the
fokker-planck equations in three dimensions which seems to quickly
exacerbate this issue.
   Thanks,
     Nathan


On Thu, May 26, 2011 at 7:31 PM, Graham Rowlands <
question159192@xxxxxxxxxxxxxxxxxxxxx> wrote:

> Question #159192 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/159192
>
> Description changed to:
> Greetings,
>
> I'm solving the Fokker-Planck equation in 2D just fine in the time-
> dependent case (Crank-Nicolson using dolfin in the usual manner I've
> seen around here.) I'm now trying to find the steady state solution
> without waiting for the dynamic simulation to converge. This is
> apparently much more difficult for me to wrap my head around.
>
> The equation (with implied summation over repeated indices):
> du/dt = d/dx_i (A_i * u) + d/dx_i (d/dx_j (B^2_ij * u ) ))
>
> maps to the following ufl code in weak form in the steady state:
> drift          = -dot(A,grad(v))*u
> diffusion  = -dot(dot(B,div(B)),grad(v))*u +
> 0.5*dot(dot(B*B.T,grad(u)),grad(v))
>
> a = (drift + diffusion)*dx  # bilinear
> L = f*v*dx                            # Linear, where f is set to 0.0 in
> the cpp code
>
> I seem to just get the trivial solution when trying the
> VariationalProblem(a,L,bcs) approach, and I guess I'm rather unsure how
> to work on relaxing this system in a similar way to that seen in the
> non-linear approach.
>
> Many thanks to all involved,
> Graham
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.