Subject: [Syfi-dev] [HG syfi] minor updates

From: [email-masked] (SyFi)
Date: Sat, 18 Apr 2009 15:24:16 +0200 (CEST)

One or more new changesets pushed to the primary SyFi repository.
A short summary of the last three changesets is included below.

changeset:   497:358dab8b5fe2245ce5da4aed35f14ee798d723ec
tag:         tip
user:        "Kent-Andre Mardal <[email-masked]>"
date:        Sat Apr 18 15:24:09 2009 +0200
files:       demo/NonLinearVectorPoisson/python/demo.py
description:
minor updates

(these test did not reveal anything other than that
 things seem to work in these simple cases)

changeset:   496:1a16250dacedbc5b4e555e629119476fae67dee0
user:        "Kent-Andre Mardal <[email-masked]>"
date:        Sat Apr 18 15:08:23 2009 +0200
files:       demo/NonLinearPoisson/python/demo.py
description:
updated demo.py

changeset:   495:032d75c80cd04bb98b347b42d856477c9140556e
user:        "Kent-Andre Mardal <[email-masked]>"
date:        Sat Apr 18 15:07:44 2009 +0200
files:       demo/NonLinearPoisson/python/manufactored_solution.py
description:
added manufactored solution script

----------------------------------------------------------------------
For more details, visit http://www.fenics.org/hg/syfi

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 15:38:06 +0200

>  things seem to work in these simple cases)

Are you seeing problems with convergence? Or does this problem converge
to an incorrect answer?

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 15:48:47 +0200 (CEST)

The nonlinear hyper-elasticity demos do not seem to produce
the correct solution. I have been looking in particular on the Fung stuff.

Newton seems to converge with quadratic convergence whenever
it is close to the solution so I think the differentiation is correct.

Kent

Kent

> _______________________________________________
> SyFi-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/syfi-dev
>

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 16:06:56 +0200

> it is close to the solution so I think the differentiation is correct.

OK. I have been trying a few things over the past week (compiling forms
with FFC instead) and here is what I see:

1. Nonlinear scalar Poisson works.
2. Linear elasticity posed as a special case of a nonlinear problem
works. (Converges in one iteration).
3. Nonlinear elasticity (hyperelasticity with the St. Venant-Kirchhoff
model) a. does not converge or b. the direct solvers complain about a
singular matrix.

I didn't think to look into the boundary conditions for 3 because I used
the same setting as 2. I will look at them more closely.

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 16:52:52 +0200 (CEST)

> the same setting as 2. I will look at them more closely.

I don't think it is a problem with the bc. I tested some nonlinear variants
of Poisson that come from convex functionals. In these cases Newton
should work. Can we say the same about SVK or Fung or do we either
need a relaxation scheme or start close to the solution ?

Kent

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 17:01:08 +0200

> need a relaxation scheme or start close to the solution ?

For the SVK or Fung (or any material model) that has been correctly
linearised, the Newton-Raphson scheme should converge quadratically as
long as one starts close to the solution. In practice, this is ensured
by not driving the problem too much in a given (time) step.

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 17:03:53 +0200 (CEST)

> by not driving the problem too much in a given (time) step.

The Newton-Raphson scheme should always converge quadratic as long as
one starts close enough to the solution :)

Kent

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 17:06:14 +0200

> one starts close enough to the solution :)

I know :). Then I didn't understand what you meant by:

I tested some nonlinear variants of Poisson that come from convex
functionals. In these cases Newton should work.

Did you mean, a Newton solver should work from any starting point?

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 17:16:29 +0200 (CEST)

> Did you mean, a Newton solver should work from any starting point?

Yes,  for an equation like
-Delta u + u**p

where p is odd.

Kent

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 17:21:51 +0200

> Yes,  for an equation like

OK, and by "work" do you mean will head in the right direction (not
necessarily at a quadratic rate)? If that is what you meant, then this
is true of these hyperelastic materials as well. The strain energy is a
convex function of the deformation measure.

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 17:32:39 +0200 (CEST)

>>>>>>>> the correct solution. I have been looking in particular on the
>>>>>>>> Fung
>>>>>>>> stuff.
>>>>>>>>

For the above equation I think there should be only one minimum and so
it will head in the right direction. I'm not sure whether it should be
quadratic.
If it is convex it can have several minium and it may therefore head
towards different minima at different times but I didn't think it could head
off to infinity (or creating a singular matrix) ?

Kent

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 17:38:06 +0200

> it will head in the right direction. I'm not sure whether it should be

That was what I was trying to clarify in an earlier e-mail. It will
always work, but converge quadratically only if close.

> off to infinity (or creating a singular matrix) ?

It should not. There is something else going on. I need to look at
things more carefully.

Harish

---

From: Jed Brown <[email-masked]>
Date: Sat, 18 Apr 2009 18:06:57 +0200

> If it is convex it can have several minium and it may therefore head
        ^^^^^^^^^

You probably mean "is not convex".

> off to infinity (or creating a singular matrix) ?

Convexity is not sufficient to guarantee global convergence, you need a
minimum to exist (consider exp(x)) and the Jacobian should be Lipshitz
continuous with J^{-1} bounded.  These conditions are probably satisfied
for the models you are discussing.

Jed

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 18:30:14 +0200

> for the models you are discussing.

These conditions are satisfied for this class of models.

Harish

---

From: [email-masked]
Date: Sat, 18 Apr 2009 19:43:30 +0200 (CEST)

> for the models you are discussing.

Right, but isn't there some results that say that if the functional F(u) goes
to infinity when u goes to +/- infinity and the functional is convex. I
probably
need to refresh this stuff.

Anyway, the main point for me was to check that the AD works properly in an
example where I know it should work. And to include non-homogenous Dirichlet
BC. It appears to work in these simple examples.

The SVK and Fung demos do still not work properly, but it might be that we
start with functions to far away. Or there might be a bug...

Kent

---

From: Jed Brown <[email-masked]>
Date: Sat, 18 Apr 2009 20:18:46 +0200

On Sat 2009-04-18 19:43, [email-masked] wrote:
> to infinity when u goes to +/- infinity and the functional is convex.

I was probably being overly pedantic, but there is something here
(although it might not apply to your problem).  If F is convex and F \to
\infty as |u| \to \infty, there is a unique minimum u^*.  This really is
not enough to for ||J^{-1}|| to be bounded.  Consider minimizing the
convex functional

  x^2, x<1
  1+x, x>1

If you start at x=2, J is singular.  J is also not Lipschitz continuous,
but mollification would make it so, boundedness of ||J^{-1}|| is still
needed.

> start with functions to far away. Or there might be a bug...

Is the Jacobian correct?  You can check this by seeing if the Krylov
iteration with matrix-free Jacobian and LU preconditioner converges in
one iteration (-snes_mf_operator -pc_type lu in PETSc).  If the Jacobian
is correct but residuals increase, you might need a line search.  I
think ||J^{-1}|| is bounded for these problems so I don't think line
search will fail.

Jed

---

From: Jed Brown <[email-masked]>
Date: Sat, 18 Apr 2009 20:23:57 +0200

> not enough to for ||J^{-1}|| to be bounded.  Consider minimizing the
    ^^^^

Should be 2x-1, but you probably got the point.

Jed

---

From: Harish Narayanan <[email-masked]>
Date: Sat, 18 Apr 2009 20:57:58 +0200

> start with functions to far away. Or there might be a bug...

I'm having problems even when starting arbitrarily close to the
solution. UMFPACK's direct solver reports that the matrix being solved
is singular.

The problem I am considering has non-homogeneous Dirichlet and
homogeneous Neumann boundary conditions. I've checked these with an
analogous linear problem and they appear to work.

Harish