Problem with variational derivative

Asked by Krishna Garikipati on 2010-06-08

I am solving a coupled problem: Nonlinear advection-diffusion and hyperelasticity

#Function spaces are:
Q = FunctionSpace(mesh, "CG", 1);
V = VectorFunctionSpace(mesh, "CG", 1);
U = VectorFunctionSpace(mesh, "CG", 1);

#Functions for advection-diff
u1c = Function(Q) #Solution from current time step
u10 = Function(Q) #Initial condition

#Functions for hyperelasticity
v2 = TestFunction(U) # Test function
du2 = TrialFunction(U) # Incremental displacement
u2c = Function(U) # Displacement in current time step

#Concentration field influencing the stress. Strain energy function psi(C), where C is the right Cauchy-Green tensor.
psi = ((u10/u1c)**1.333333)*const1*(tr(C))**2 - ((u10/u1c)**0.666666)*const2*tr(C) + ((u1c/u10)**1.333333)*const3*tr(C*C)
P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) #Stress

#Variational form
L2 = inner(P, grad(v2))*dx(1) - inner(B, v2)*dx(1) - inner(T2, v2)*ds(1)
a2 = derivative(L2, u2c, du2)
# Solve nonlinear variational problem
problem = VariationalProblem(a2, L2, [bclm, bcrm],
                             cell_domains=sub_domains,
                             exterior_facet_domains=boundary,
                             nonlinear=True)

This code takes inordinately long to form the matrices, but appears to converge quadratically to the correct solution. If instead I move the multiplication by (u10/u1c) from
P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) to the previous line

psi = ((u10/u1c))*const1*(tr(C))**2 - ((u10/u1c)**0.333333)*const2*tr(C) + ((u1c/u10))*const3*tr(C*C)
P = 2*F*diff(psi, C) #Stress

the matrices are assembled faster by a factor of 30 or 40, but the starting residual is different and I lose quadratic convergence of the residual.

Is there something obviously wrong? It is unexpected that rewriting the stress without changing its mathematical form causes such a difference in the assembly and solution procedures.

 ----Krishna

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2010-06-08
Last reply:
2010-06-08
Garth Wells (garth-wells) said : #1

On 08/06/10 16:19, Krishna Garikipati wrote:
> New question #113926 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/113926
>
> I am solving a coupled problem: Nonlinear advection-diffusion and hyperelasticity
>
> #Function spaces are:
> Q = FunctionSpace(mesh, "CG", 1);
> V = VectorFunctionSpace(mesh, "CG", 1);
> U = VectorFunctionSpace(mesh, "CG", 1);
>
> #Functions for advection-diff
> u1c = Function(Q) #Solution from current time step
> u10 = Function(Q) #Initial condition
>
> #Functions for hyperelasticity
> v2 = TestFunction(U) # Test function
> du2 = TrialFunction(U) # Incremental displacement
> u2c = Function(U) # Displacement in current time step
>
> #Concentration field influencing the stress. Strain energy function psi(C), where C is the right Cauchy-Green tensor.
> psi = ((u10/u1c)**1.333333)*const1*(tr(C))**2 - ((u10/u1c)**0.666666)*const2*tr(C) + ((u1c/u10)**1.333333)*const3*tr(C*C)
> P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) #Stress
>
> #Variational form
> L2 = inner(P, grad(v2))*dx(1) - inner(B, v2)*dx(1) - inner(T2, v2)*ds(1)
> a2 = derivative(L2, u2c, du2)
> # Solve nonlinear variational problem
> problem = VariationalProblem(a2, L2, [bclm, bcrm],
> cell_domains=sub_domains,
> exterior_facet_domains=boundary,
> nonlinear=True)
>
>
> This code takes inordinately long to form the matrices, but appears to converge quadratically to the correct solution.

Try putting

     parameters["form_compiler"]["cpp_optimize"] = True
     parameters["form_compiler"]["optimize"] = True

somewhere near the top of your file. It will take longer to compile the
code the first time through, but the assembly should be much faster
(likely orders of magnitude faster).

> If instead I move the multiplication by (u10/u1c) from
> P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) to the previous line
>
> psi = ((u10/u1c))*const1*(tr(C))**2 - ((u10/u1c)**0.333333)*const2*tr(C) + ((u1c/u10))*const3*tr(C*C)
> P = 2*F*diff(psi, C) #Stress
>
> psi = ((u10/u1c)**1.333333)*const1*(tr(C))**2 - ((u10/u1c)**0.666666)*const2*tr(C) + ((u1c/u10)**1.333333)*const3*tr(C*C)
> P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) #Stress

> the matrices are assembled faster by a factor of 30 or 40, but the starting residual is different and I lose quadratic convergence of the residual.
>
> Is there something obviously wrong? It is unexpected that rewriting the stress without changing its mathematical form causes such a difference in the assembly and solution procedures.
>

It should of course give the same result if the expressions are the
same, but if may affect the runtime, particularly when optimisations are
switched off (Kristian could comment in more details).

Garth

> ----Krishna
>
>
>

On 8 June 2010 17:19, Krishna Garikipati
<email address hidden> wrote:
> New question #113926 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/113926
>
> I am solving a coupled problem: Nonlinear advection-diffusion and hyperelasticity
>
> #Function spaces are:
> Q = FunctionSpace(mesh, "CG", 1);
> V = VectorFunctionSpace(mesh, "CG", 1);
> U = VectorFunctionSpace(mesh, "CG", 1);
>
> #Functions for advection-diff
> u1c  = Function(Q)          #Solution from current time step
> u10  = Function(Q)          #Initial condition
>
> #Functions for hyperelasticity
> v2   = TestFunction(U)       # Test function
> du2  = TrialFunction(U)      # Incremental displacement
> u2c  = Function(U)           # Displacement in current time step
>
> #Concentration field influencing the stress. Strain energy function psi(C), where C is the right Cauchy-Green tensor.
> psi = ((u10/u1c)**1.333333)*const1*(tr(C))**2 - ((u10/u1c)**0.666666)*const2*tr(C) + ((u1c/u10)**1.333333)*const3*tr(C*C)
> P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) #Stress
>
> #Variational form
> L2 = inner(P, grad(v2))*dx(1) - inner(B, v2)*dx(1) - inner(T2, v2)*ds(1)
> a2 = derivative(L2, u2c, du2)

Could you send the entire definition of the variational form? Some
variables like const1, const2, const3, B, T2, C (did I miss any?) are
not defined, I can't get FFC to compile without.

> # Solve nonlinear variational problem
> problem = VariationalProblem(a2, L2, [bclm, bcrm],
>                             cell_domains=sub_domains,
>                             exterior_facet_domains=boundary,
>                             nonlinear=True)
>
>
> This code takes inordinately long to form the matrices, but appears to converge quadratically to the correct solution. If instead I move the multiplication by (u10/u1c) from
> P = 2*((u1c/u10)**0.333333)*F*diff(psi, C) to the previous line
>
> psi = ((u10/u1c))*const1*(tr(C))**2 - ((u10/u1c)**0.333333)*const2*tr(C) + ((u1c/u10))*const3*tr(C*C)
> P = 2*F*diff(psi, C) #Stress
>
> the matrices are assembled faster by a factor of 30 or 40, but the starting residual is different and I lose quadratic convergence of the residual.

The expressions for psi and P have been simplified since only one
power (**0.3333) remains. These kind of reductions are not done in UFL
or FFC so that might explain the difference in runtime.
The results should of course be the same, but as they are not this can
of course be another reason why the runtime differs (because you
assemble something different).
To get to the bottom of this I guess we need to compare the UFL
expressions or look at the generated code.

Kristian

> Is there something obviously wrong? It is unexpected that rewriting the stress without changing its mathematical form causes such a difference in the assembly and solution procedures.
>
>  ----Krishna
>
>
>
> --
> 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     : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>

Can you help with this problem?

Provide an answer of your own, or ask Krishna Garikipati for more information if necessary.

To post a message you must log in.