StateMaps

struct PDEConstrainedFunctionals{N,A} <: AbstractPDEConstrainedFunctionals{N}

An object that computes the objective, constraints, and their derivatives.

Implementation

This implementation computes derivatives of a integral quantity

$F(u(\varphi),\varphi,\mathrm{d}\Omega_1,\mathrm{d}\Omega_2,...) = \Sigma_{i}\int_{\Omega_i} f_i(\varphi)~\mathrm{d}\Omega$

with respect to an auxiliary parameter $\varphi$ where $u$ is the solution to a PDE and implicitly depends on $\varphi$. This requires two pieces of information:

1. Computation of $\frac{\partial F}{\partial u}$ and $\frac{\partial F}{\partial \varphi}$ (handled by StateParamMap).

2. Computation of $\frac{\partial F}{\partial u} \frac{\partial u}{\partial \varphi}$ at $\varphi$ and $u$ using the adjoint method (handled by AbstractFEStateMap). I.e., let

   $\frac{\partial F}{\partial u} \frac{\partial u}{\partial \varphi} = -\lambda^\intercal \frac{\partial \mathcal{R}}{\partial \varphi}$

   where $\mathcal{R}$ is the residual and solve the (linear) adjoint problem:

   $\frac{\partial \mathcal{R}}{\partial u}^\intercal\lambda = \frac{\partial F}{\partial u}^\intercal.$

The gradient is then $\frac{\partial F}{\partial \varphi} = \frac{\partial F}{\partial \varphi} - \frac{\partial F}{\partial u}\frac{\partial u}{\partial \varphi}$.

Parameters

• J: A StateParamMap corresponding to the objective.
• C: A vector of StateParamMap corresponding to the constraints.
• dJ: The DoFs for the objective sensitivity.
• dC: The DoFs for each constraint sensitivity.
• analytic_dJ: a Function for computing the analytic objective sensitivity.
• analytic_dC: A vector of Function for computing the analytic objective sensitivities.
• state_map::A: The state map for the problem.

Note

• If analytic_dJ = nothing automatic differentiation will be used.
• If analytic_dC[i] = nothing automatic differentiation will be used for C[i].

Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh)

Evaluate the objective and constraints, and their derivatives at φh.

evaluate_functionals!(pcf::PDEConstrainedFunctionals,φh)

Evaluate the objective and constraints at φh.

evaluate_derivatives!(pcf::PDEConstrainedFunctionals,φh)

Evaluate the derivatives of the objective and constraints at φh.

get_state(m::AbstractFEStateMap)

Return the solution/state u to the FE problem.

struct StateParamMap{A,B,C,D} <: AbstractStateParamMap

A wrapper to handle partial differentation of a function F of a specific form (see below) in a ChainRules.jl compatible way with caching.

Assumptions

We assume that we have a function F of the following form:

F(u,φ) = ∫(f(u,φ))dΩ₁ + ∫(g(u,φ))dΩ₂ + ...,.

where u and φ are each expected to inherit from Union{FEFunction,MultiFieldFEFunction} or the GridapDistributed equivalent.

ChainRulesCore.rrule(u_to_j::StateParamMap,uh,φh)

Return the evaluation of a StateParamMap and a a function for evaluating the pullback of u_to_j. This enables compatiblity with ChainRules.jl

abstract type AbstractFEStateMap

Types inheriting from this abstract type should enable the evaluation and differentiation of the solution to an FE problem u that implicitly depends on an auxiliary parameter φ.

struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap

A structure to enable the forward problem and pullback for affine finite element operators AffineFEOperator.

Parameters

• biform::A: Function defining the bilinear form.
• liform::B: Function defining the linear form.
• spaces::C: Tuple of finite element spaces.
• plb_caches::D: A cache for the pullback operator.
• fwd_caches::E: A cache for the forward problem.
• adj_caches::F: A cache for the adjoint problem.

AffineFEStateMap(
  a::Function,l::Function,
  U,V,V_φ,φh;
  assem_U = SparseMatrixAssembler(U,V),
  assem_adjoint = SparseMatrixAssembler(V,U),
  assem_deriv = SparseMatrixAssembler(V_φ,V_φ),
  ls::LinearSolver = LUSolver(),
  adjoint_ls::LinearSolver = LUSolver()
)

Create an instance of AffineFEStateMap given the bilinear form a and linear form l as Function types, trial and test spaces U and V, the FE space V_φ for φh and derivatives, and the measures as additional arguments.

Optional arguments enable specification of assemblers and linear solvers.

struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap

A structure to enable the forward problem and pullback for nonlinear finite element operators.

Parameters

• res::A: a Function defining the residual of the problem.
• jac::B: a Function defining Jacobian of the residual.
• spaces::C: Tuple of finite element spaces.
• plb_caches::D: A cache for the pullback operator.
• fwd_caches::E: A cache for the forward problem.
• adj_caches::F: A cache for the adjoint problem.

NonlinearFEStateMap(
  res::Function,jac::Function,adjoint_jac::Function,U,V,V_φ,φh;
  assem_U = SparseMatrixAssembler(U,V),
  assem_adjoint = SparseMatrixAssembler(V,U),
  assem_deriv = SparseMatrixAssembler(V_φ,V_φ),
  nls::NonlinearSolver = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-8,verbose=true),
  adjoint_ls::LinearSolver = LUSolver()
)

In addition to the above, pass the jacobian adjoint_jac for the purposes of solving the adjoint problem. This can be computed with AD or by hand, but and allows the user to specify a different jacobian for the forward problem (e.g., for picard iterations).

NonlinearFEStateMap(
  res::Function,jac::Function,U,V,V_φ,φh;
  assem_U = SparseMatrixAssembler(U,V),
  assem_adjoint = SparseMatrixAssembler(V,U),
  assem_deriv = SparseMatrixAssembler(V_φ,V_φ),
  nls::NonlinearSolver = NewtonSolver(LUSolver();maxiter=50,rtol=1.e-8,verbose=true),
  adjoint_ls::LinearSolver = LUSolver()
)

Create an instance of NonlinearFEStateMap given the residual res as a Function type, trial and test spaces U and V, the FE space V_φ for φh and derivatives, and the measures as additional arguments.

Optional arguments enable specification of assemblers, nonlinear solver, and adjoint (linear) solver.

struct RepeatingAffineFEStateMap <: AbstractFEStateMap

A structure to enable the forward problem and pullback for affine finite element operators AffineFEOperator with multiple linear forms but only a single bilinear form.

Parameters

• biform: A Function defining the bilinear form.
• liform: A vector of Function defining the linear forms.
• spaces: Repeated finite element spaces.
• spaces_0: Original finite element spaces that are being repeated.
• plb_caches: A cache for the pullback operator.
• fwd_caches: A cache for the forward problem.
• adj_caches: A cache for the adjoint problem.

RepeatingAffineFEStateMap(
  nblocks::Int,biform::Function,liforms::Vector{<:Function},
  U0,V0,V_φ,φh;
  assem_U = SparseMatrixAssembler(U0,V0),
  assem_adjoint = SparseMatrixAssembler(V0,U0),
  assem_deriv = SparseMatrixAssembler(V_φ,V_φ),
  ls::LinearSolver = LUSolver(),
  adjoint_ls::LinearSolver = LUSolver()
)

Create an instance of RepeatingAffineFEStateMap given the number of blocks nblocks, a bilinear form a, a vector of linear form l as Function types, the trial and test spaces U and V, the FE space V_φ for φh and derivatives, and the measures as additional arguments.

Optional arguments enable specification of assemblers and linear solvers.

Note

• The resulting FEFunction will be a MultiFieldFEFunction (or GridapDistributed equivalent) where each field corresponds to an entry in the vector of linear forms

struct StaggeredAffineFEStateMap{NB,SB} <: AbstractFEStateMap{NB,SB}
  biforms    :: Vector{<:Function}
  liforms    :: Vector{<:Function}
  ∂Rk∂xhi    :: Tuple{Vararg{Tuple{Vararg{Function}}}}
  spaces     :: A
  assems     :: B
  solvers    :: C
  plb_caches :: D
  fwd_caches :: E
  adj_caches :: F
end

Affine staggered state map for the equivalent StaggeredAffineFEOperator, used to solve staggered problems where the k-th equation is linear in u_k.

Similar to the StaggeredAffineFEOperator counterpart, we expect a set of bilinear/linear form pairs that also depend on φ:

a_k((u_1,...,u_{k-1}),u_k,v_k,φ) = ∫(...)
l_k((u_1,...,u_{k-1}),v_k,φ) = ∫(...)

These can be assembled into a set of linear systems:

A_k u_k = b_k

where A_k and b_k only depend on the previous variables u_1,...,u_{k-1}.

## Weak forms
a1(u1,v1,φ) = ...
l1(v1,φ) = ...
# Treat (u1,φ) as the primal variable
a2(u2,v2,(u1,φ)) = ...
l2(v2,(u1,φ)) = ...

## Build StateMaps
φ_to_u1 = AffineFEStateMap(a1,l1,U1,V,V_φ,φh)
# u1φ_to_u2 has a MultiFieldFESpace V_u1φ of primal vars
u1φ_to_u2 = AffineFEStateMap(a2,l2,U2,V,V_u1φ,interpolate([1,φh],V_u1φ))
# The StateParamMap F needs to take a MultiFieldFEFunction u1u2h ∈ U_u1u2
F = GridapTopOpt.StateParamMap(F,U_u1u2,V_φ,assem_U_u1u2,assem_V_φ)

function φ_to_j(φ)
  u1 = φ_to_u1(φ)
  u1φ = combine_fields(V_u1φ,u1,φ) # Combine vectors of DOFs
  u2 = u1φ_to_u2(u1φ)
  u1u2 = combine_fields(U_u1u2,u1,u2)
  F(u1u2,φ)
end

pcf = CustomPDEConstrainedFunctionals(...)

StaggeredAffineFEStateMap(
    op              :: StaggeredAffineFEOperator{NB,SB},
    V_φ,
    φh;
    assem_deriv     :: Assembler = SparseMatrixAssembler(V_φ,V_φ),
    assems_adjoint  :: Vector{<:Assembler} = map(SparseMatrixAssembler,op.tests,op.trials),
    solver          :: StaggeredFESolver{NB} = StaggeredFESolver(fill(LUSolver(),length(op.biforms))),
    adjoint_solver  :: StaggeredFESolver{NB} = StaggeredFESolver(fill(LUSolver(),length(op.biforms)))
) where {NB,SB}

Create an instance of StaggeredAffineFEStateMap given a StaggeredAffineFEOperator op, the auxiliary space V_φ for φh and derivatives, and the parameter φh.

Otional arguemnts:

• assem_deriv is the assembler for the derivative space.
• assems_adjoint is a vector of assemblers for the adjoint space.
• solver is a StaggeredFESolver for the forward problem.
• adjoint_solver is a StaggeredFESolver for the adjoint problem.

mutable struct StaggeredNonlinearFEStateMap{NB,SB} <: AbstractFEStateMap{NB,SB}
  const residuals         :: Vector{<:Function}
  const jacobians         :: Vector{<:Function}
  const adjoint_jacobians :: Vector{<:Function}
  const ∂Rk∂xhi           :: Tuple{Vararg{Tuple{Vararg{Function}}}}
  const spaces            :: A
  const assems            :: B
  const solvers           :: C
  const plb_caches        :: D
  fwd_caches              :: E
  const adj_caches        :: F
end

Staggered nonlinear state map for the equivalent StaggeredNonlinearFEOperator, used to solve staggered problems where the k-th equation is nonlinear in u_k.

Similar to the previous structure and the StaggeredNonlinearFEOperator counterpart, we expect a set of residual/jacobian pairs that also depend on φ:

jack((u1,...,u{k-1},φ),uk,duk,dvk) = ∫(...) resk((u1,...,u{k-1},φ),uk,v_k) = ∫(...)

function StaggeredNonlinearFEStateMap(
  op                :: StaggeredNonlinearFEOperator{NB,SB},
  V_φ,
  φh;
  assem_deriv       :: Assembler = SparseMatrixAssembler(V_φ,V_φ),
  assems_adjoint    :: Vector{<:Assembler} = map(SparseMatrixAssembler,op.tests,op.trials),
  solver            :: StaggeredFESolver{NB} = StaggeredFESolver(
    fill(NewtonSolver(LUSolver();maxiter=50,rtol=1.e-8,verbose=true),length(op.residuals))),
  adjoint_solver    :: StaggeredFESolver{NB} = StaggeredFESolver(fill(LUSolver(),length(op.residuals))),
  adjoint_jacobians :: Vector{<:Function} = op.jacobians
) where {NB,SB}

Create an instance of StaggeredNonlinearFEStateMap given a StaggeredNonlinearFEOperator op, the auxiliary space V_φ for φh and derivatives, and the parameter φh.

Otional arguemnts:

• assem_deriv is the assembler for the derivative space.
• assems_adjoint is a vector of assemblers for the adjoint space.
• solver is a StaggeredFESolver for the forward problem.
• adjoint_solver is a StaggeredFESolver for the adjoint problem.
• adjoint_jacobians is a vector of jacobians for the adjoint problem.

rrule(φ_to_u::AbstractFEStateMap,φh)

Return the evaluation of a AbstractFEStateMap and a a function for evaluating the pullback of φ_to_u. This enables compatiblity with ChainRules.jl

pullback(φ_to_u::AbstractFEStateMap,uh,φh,du;updated)

Compute ∂F∂u*dudφ at φh and uh using the adjoint method. I.e., let

∂F∂u*dudφ = -λᵀ*dRdφ

and solve the adjoint problem dRduᵀ*λ = ∂F∂uᵀ using adjoint_solve!.

forward_solve!(φ_to_u::AbstractFEStateMap,φh)

Evaluate the forward problem u given φ. This should compute the FE problem.

adjoint_solve!(φ_to_u::AbstractFEStateMap,du::AbstractVector)

Evaluate the solution to the adjoint problem given a RHS vector ∂F∂u denoted du. This should solve the linear problem dRduᵀ*λ = ∂F∂uᵀ.

update_adjoint_caches!(φ_to_u::AbstractFEStateMap,uh,φh)

Update the cache for the adjoint problem. This is usually a tuple of objects.

dRdφ(φ_to_u::AbstractFEStateMap,uh,vh,φh)

Compute the derivative with respect to φh of the residual R.

get_state(m::AbstractFEStateMap)

Return the solution/state u to the FE problem.

get_spaces(m::AbstractFEStateMap)

Return a collection of FE spaces. The first four entires should correspond to get_trial_space, get_test_space, get_aux_space, and get_deriv_space unless these are overloaded for a particular implementation.

get_assemblers(m::AbstractFEStateMap)

Return a collection of assemblers. The first two entires should correspond to get_pde_assembler and get_deriv_assembler unless these are overloaded for a particular implementation.

get_trial_space(m::AbstractFEStateMap)

Return trial space for FE problem.

get_test_space(m::AbstractFEStateMap)

Return test space for FE problem.

get_aux_space(m::AbstractFEStateMap)

Return space for auxillary parameter.

get_deriv_space(m::AbstractFEStateMap)

Return space for derivatives.

get_pde_assembler(m::AbstractFEStateMap)

Return assembler for FE problem.

get_deriv_assembler(m::AbstractFEStateMap)

Return assembler for derivatives.

Gridap.jacobian(F,uh::Vector{<:CellField},K::Int)

Given a function F that returns a DomainContribution when called, and a vector of FEFunctions or CellField uh, evaluate the Jacobian F with respect to uh[K].

Gridap.gradient(F,uh::Vector{<:CellField},K::Int)

Given a function F that returns a DomainContribution when called, and a vector of FEFunctions uh, evaluate the partial derivative of F with respect to uh[K].

Example

Suppose uh and φh are FEFunctions with measures dΩ and dΓ_N. Then the partial derivative of a function J wrt to φh is computed via

J(u,φ) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N
∂J∂φh = ∇(J,[uh,φh],2)

where f and g are user defined.