dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

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