parameter dependant model

[rveltz]

Hi,

This is not really an issue. I would like to hear your view on making the pde formulations parameter dependent. I am using your package you my bifurcation tools. I have written a simple wrapper to use my "newton".

I guess the model parameters could be put in your cache from _new_nlsolve_cache.

I would like to hear your opinion about this and potential pitfalls.

Thank you,

[fverdugo]

Hi @rveltz ,

nice that you have interested in the project.

I don't really understand what you mean by "parameter dependent". Are you willing to include state variables in your PDE?

This should be possible in Gridap (even tough it is a quite experimental feature at the moment). See for instance the isotropic damage tutorial.

[rveltz]

It is somehow like the example you mention but I would like to avoid using a closure or const variables.

By the way, the example you show, the load is on each small side? Is it an Euler buckling rod?

Say you have the PDE 0=\Delta u + \lambda (e^u-u) with Neumann BC. I want to solve this PDE for different values of \lambda.

Note that I want to use PseudoArcLengthContinuation.jl, so I will use its newton solver and continuation. For using continuation, I need to specify a parameter dependent PDE in the form (out of place) F(u, p) where p is the set of parameters. This is easy to do once you have access to your cache from _new_nlsolve_cache. Hence the question is about recomputing this cache for each parameter value I guess.

In the case of the PDE above, calling _new_nlsolve_cache for each lambda would be inefficient because the discretization basically does not change (same BC, same jacobian sparsity pattern...).

[fverdugo]

@rveltz As you say, in the example I have provided, we recompute from scratch everything (e.g. the jacopian). However, this can be easily improved by properly using the API of our solvers.

uh1, cache  = solve!(uh0, solver, op1)
uh2, cache  = solve!(uh1, solver, op2, cache) # Reuse cache from previous run

This works as long as the sparsity pattern of the operators op1 and op2 is the same. Each specialization of solver defines what is stored in the cache. The NLSolver specialization reuses the sparsity pattern of the jacobian among other things.

By the way, the example you show, the load is on each small side? Is it an Euler buckling rod?

In this case, we are using 3D continuum elements to model the beam, i.e., small displacement linear elasticity, with a non-linear state-dependent constitutive law. That is, we don't have buckling, the non-linearity only comes from the material law. The small sides are clamped and the beam bends under self weight.

By the way, which quantities do you need for your arc length conitnuation? F(u,p) and dF_du(u,p) (the jacobian wrt u) can be computed with Gridap.

[fverdugo]

Provably you also need dF_dp(u,p) (jacobian wrt p). I think computing this is also possible.

In any case, I believe that all what you need can be done now, but it requires some advanced expertise with the library.

[rveltz]

This is what I do for newton

function _new_nlsolve_cache_RV(u,nls,feop)
	@assert Gridap.FESpaces.is_a_fe_function(u)
	x0 = get_free_values(u)
	op = Gridap.FESpaces.get_algebraic_operator(feop)
	f!(r,x) = Gridap.FESpaces.residual!(r,op,x)
	j!(j,x) = Gridap.FESpaces.jacobian!(j,op,x)
	fj!(r,j,x) = residual_and_jacobian!(r,j,op,x)
	f0, j0 = Gridap.FESpaces.residual_and_jacobian(op,x0)
	ss = symbolic_setup(nls.ls,j0)
	ns = numerical_setup(ss,j0)
	return (x0 = x0, f0 = f0, j0 = j0, f = f!, df = j!, ns = ns)
end
	_new_nlsolve_cache_RV(uh0, nls, op)

# I pass the cache as a parameter for now
function Fbratu(x, par)
	out = similar(x)
	par.f(out, x)
	return out
end

function Jbratu(x, par)
	par.df(par.j0, x)
	return par.j0
end

# linear solver to use PseudoArcLengthContinuation.jl
struct GridapLinSolve{T} <: PALC.AbstractLinearSolver
	ns::T
end

# we could save x in a struct to avoid allocations
# let's keep it simple for now
function (ls::GridapLinSolve)(A, rhs)
	x = similar(rhs)
	numerical_setup!(ls.ns,A)
	solve!(x,ls.ns,rhs)
	return x, true, 1
end

# cache for solvers
cache_palc = _new_nlsolve_cache_RV((uh0), nls, op)
# linear solver to use "my" newton
ls  = GridapLinSolve(cache_palc.ns)
# newton options
optnew = NewtonPar(verbose = true, linsolver = ls, maxIter = 10)
# Newton Krylov solver
PALC.newton(Fbratu, Jbratu, cache_palc.x0, cache_palc, optnew, norm = x->norm(x,Inf))

By the way, which quantities do you need for your arc length conitnuation? F(u,p) and dF_du(u,p) (the jacobian wrt u) can be computed with Gridap.

Yes, this is what you need. See above for an example where the parameters are actually not used.

Provably you also need dF_dp(u,p) (jacobian wrt p). I think computing this is also possible.

This is computed by finite differences for now, it does not need to be that precise in practice.

[fverdugo]

Hi @rveltz and @santiagobadia

Here, I provide a small example on how to eficiently solve a parameter-dependent problem by reusing as much data as possible between the solution of different parameters.

using Gridap

function main(;n)

  domain = (0,1,0,1)
  cells = (n,n)
  model = CartesianDiscreteModel(domain,cells)
  
  order = 1
  V = TestFESpace(
    model=model,reffe=:Lagrangian,valuetype=Float64,
    order=order,conformity=:H1,dirichlet_tags="boundary")
  U = TrialFESpace(V)
  
  trian = Triangulation(model)
  degree = 2*order
  quad = CellQuadrature(trian,degree)

  # Given a parameter, return a FEOperator
  function op_from_param(p)
    @law NL(u) = exp(u)
    res(u,v) = ∇(v)⋅∇(u) - p.λ * v ⋅ (NL(u) - u)
    jac(u,du,v) = ∇(v)⋅∇(du) -   p.λ * v ⋅ du ⋅ (NL(u)-1)
    t_Ω = FETerm(res,jac,trian,quad)
    op = FEOperator(U,V,t_Ω)
    op
  end

  # Setup solvers
  ls = LUSolver()
  nls = NLSolver(ls, show_trace=true, method=:newton)
  solver = FESolver(nls)

  # Initial guess and initial cache
  cache = nothing
  uh = zero(U)

  # Solve the problem for several parameters by reusing cache.
  # In particular, here we only allocate the residual and jacobian once,
  # we reuse symbolic setup of the solver (if any) etc.
  for (i,p) in enumerate([ (λ=λ,) for λ in range(0,10,length=20)])
    op = op_from_param(p)
    uh, cache = solve!(uh,solver,op,cache)
    writevtk(trian,"results_$i",cellfields=["uh"=>uh])
  end

end

main(n=10)

As you can see, this is quite straight-forward. I believe that this is the starting point in order to do the coupling with a bifurcation solver.

[fverdugo]

@rveltz, now, you can define the functions you need for bifurcation as follows:

using Gridap.FESpaces

function F(u,p)
  op = op_from_param(p)
  algop = get_algebraic_operator(op)
  r = residual(algop,u)
  r
end

function J_u(u,p)
  op = op_from_param(p)
  algop = get_algebraic_operator(op)
  j = jacobian(algop,u)
  j
end

This can be improved if your library accepts "implace" versions of these functions. Then, you would do something like:

function F!(r,u,p)
  op = op_from_param(p)
  algop = get_algebraic_operator(op)
  residual!(r,algop,u)
  r
end

function J_u!(j,u,p)
  op = op_from_param(p)
  algop = get_algebraic_operator(op)
  jacobian!(j,algop,u)
end

Computing the jacobian w.r.t p would require some more lines of code, but you say that you don't need it for the moment. If your library acceps a "matrix-free" jacobian for p, this can be also done. You would only need to create another function op_from_dparam(p,dp) in which you provide the "direction" dp in the parameter space.

[rveltz]

This is beautiful! I will play with it