# Copyright (c) 2016-2018, The University of Texas at Austin
# & University of California--Merced.
# Copyright (c) 2019-2020, The University of Texas at Austin
# University of California--Merced, Washington University in St. Louis.
#
# All Rights reserved.
# See file COPYRIGHT for details.
#
# This file is part of the hIPPYlib library. For more information and source code
# availability see https://hippylib.github.io.
#
# hIPPYlib is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License (as published by the Free
# Software Foundation) version 2.0 dated June 1991.
import dolfin as dl
import ufl
from .variables import STATE, PARAMETER, ADJOINT
from ..algorithms.linalg import Transpose
from ..algorithms.linSolvers import PETScLUSolver
from ..utils.vector2function import vector2Function
[docs]class PDEProblem(object):
""" Consider the PDE problem:
Given :math:`m`, find :math:`u` such that
.. math:: F(u, m, p) = ( f(u, m), p) = 0, \\quad \\forall p.
Here :math:`F` is linear in :math:`p`, but it may be non linear in :math:`u` and :math:`m`.
"""
[docs] def generate_state(self):
""" Return a vector in the shape of the state. """
raise NotImplementedError("Child class should implement method generate_state")
[docs] def generate_parameter(self):
""" Return a vector in the shape of the parameter. """
raise NotImplementedError("Child class should implement method generate_parameter")
[docs] def init_parameter(self, m):
""" Initialize the parameter. """
raise NotImplementedError("Child class should implement method init_parameter")
[docs] def solveFwd(self, state, x):
""" Solve the possibly nonlinear forward problem:
Given :math:`m`, find :math:`u` such that
.. math:: \\delta_p F(u, m, p;\\hat{p}) = 0, \\quad \\forall \\hat{p}.
"""
raise NotImplementedError("Child class should implement method solveFwd")
[docs] def solveAdj(self, adj, x, adj_rhs):
""" Solve the linear adjoint problem:
Given :math:`m`, :math:`u`; find :math:`p` such that
.. math:: \\delta_u F(u, m, p;\\hat{u}) = 0, \\quad \\forall \\hat{u}.
"""
raise NotImplementedError("Child class should implement method solveAdj")
[docs] def evalGradientParameter(self, x, out):
"""Given :math:`u, m, p`; evaluate :math:`\\delta_m F(u, m, p; \\hat{m}),\\, \\forall \\hat{m}.` """
raise NotImplementedError("Child class should implement method evalGradientParameter")
[docs] def setLinearizationPoint(self,x, gauss_newton_approx):
""" Set the values of the state and parameter
for the incremental forward and adjoint solvers.
Set whether Gauss Newton approximation of
the Hessian should be used."""
raise NotImplementedError("Child class should implement method setLinearizationPoint")
[docs] def solveIncremental(self, out, rhs, is_adj):
""" If :code:`is_adj = False`:
Solve the forward incremental system:
Given :math:`u, m`, find :math:`\\tilde{u}` such that
.. math::
\\delta_{pu} F(u, m, p; \\hat{p}, \\tilde{u}) = \\mbox{rhs}, \\quad \\forall \\hat{p}.
If :code:`is_adj = True`:
Solve the adjoint incremental system:
Given :math:`u, m`, find :math:`\\tilde{p}` such that
.. math::
\\delta_{up} F(u, m, p; \\hat{u}, \\tilde{p}) = \\mbox{rhs}, \\quad \\forall \\hat{u}.
"""
raise NotImplementedError("Child class should implement method solveIncremental")
[docs] def apply_ij(self,i,j, dir, out):
"""
Given :math:`u, m, p`; compute
:math:`\\delta_{ij} F(u, m, p; \\hat{i}, \\tilde{j})` in the direction :math:`\\tilde{j} =` :code:`dir`,
:math:`\\forall \\hat{i}.`
"""
raise NotImplementedError("Child class should implement method apply_ij")
[docs] def apply_ijk(self,i,j,k, x, jdir, kdir, out):
"""
Given :code:`x = [u,a,p]`; compute
:math:`\\delta_{ijk} F(u,a,p; \\hat{i}, \\tilde{j}, \\tilde{k})`
in the direction :math:`(\\tilde{j},\\tilde{k}) = (`:code:`jdir,kdir`), :math:`\\forall \\hat{i}.`
"""
raise NotImplementedError("Child class should implement apply_ijk")
[docs]class PDEVariationalProblem(PDEProblem):
def __init__(self, Vh, varf_handler, bc, bc0, is_fwd_linear = False):
self.Vh = Vh
self.varf_handler = varf_handler
if type(bc) is dl.DirichletBC:
self.bc = [bc]
else:
self.bc = bc
if type(bc0) is dl.DirichletBC:
self.bc0 = [bc0]
else:
self.bc0 = bc0
self.A = None
self.At = None
self.C = None
self.Wmu = None
self.Wmm = None
self.Wuu = None
self.solver = None
self.solver_fwd_inc = None
self.solver_adj_inc = None
self.is_fwd_linear = is_fwd_linear
self.n_calls = {"forward": 0,
"adjoint":0 ,
"incremental_forward":0,
"incremental_adjoint":0}
[docs] def generate_state(self):
""" Return a vector in the shape of the state. """
return dl.Function(self.Vh[STATE]).vector()
[docs] def generate_parameter(self):
""" Return a vector in the shape of the parameter. """
return dl.Function(self.Vh[PARAMETER]).vector()
[docs] def init_parameter(self, m):
""" Initialize the parameter. """
dummy = self.generate_parameter()
m.init( dummy.mpi_comm(), dummy.local_range() )
[docs] def solveFwd(self, state, x):
""" Solve the possibly nonlinear forward problem:
Given :math:`m`, find :math:`u` such that
.. math:: \\delta_p F(u, m, p;\\hat{p}) = 0,\\quad \\forall \\hat{p}."""
self.n_calls["forward"] += 1
if self.solver is None:
self.solver = self._createLUSolver()
if self.is_fwd_linear:
u = dl.TrialFunction(self.Vh[STATE])
m = vector2Function(x[PARAMETER], self.Vh[PARAMETER])
p = dl.TestFunction(self.Vh[ADJOINT])
res_form = self.varf_handler(u, m, p)
A_form = ufl.lhs(res_form)
b_form = ufl.rhs(res_form)
A, b = dl.assemble_system(A_form, b_form, bcs=self.bc)
self.solver.set_operator(A)
self.solver.solve(state, b)
else:
u = vector2Function(x[STATE], self.Vh[STATE])
m = vector2Function(x[PARAMETER], self.Vh[PARAMETER])
p = dl.TestFunction(self.Vh[ADJOINT])
res_form = self.varf_handler(u, m, p)
dl.solve(res_form == 0, u, self.bc)
state.zero()
state.axpy(1., u.vector())
[docs] def solveAdj(self, adj, x, adj_rhs):
""" Solve the linear adjoint problem:
Given :math:`m, u`; find :math:`p` such that
.. math:: \\delta_u F(u, m, p;\\hat{u}) = 0, \\quad \\forall \\hat{u}.
"""
self.n_calls["adjoint"] += 1
if self.solver is None:
self.solver = self._createLUSolver()
u = vector2Function(x[STATE], self.Vh[STATE])
m = vector2Function(x[PARAMETER], self.Vh[PARAMETER])
p = dl.Function(self.Vh[ADJOINT])
du = dl.TestFunction(self.Vh[STATE])
dp = dl.TrialFunction(self.Vh[ADJOINT])
varf = self.varf_handler(u, m, p)
adj_form = dl.derivative( dl.derivative(varf, u, du), p, dp )
Aadj, dummy = dl.assemble_system(adj_form, ufl.inner(u,du)*ufl.dx, self.bc0)
self.solver.set_operator(Aadj)
self.solver.solve(adj, adj_rhs)
[docs] def evalGradientParameter(self, x, out):
"""Given :math:`u, m, p`; evaluate :math:`\\delta_m F(u, m, p; \\hat{m}),\\, \\forall \\hat{m}.` """
u = vector2Function(x[STATE], self.Vh[STATE])
m = vector2Function(x[PARAMETER], self.Vh[PARAMETER])
p = vector2Function(x[ADJOINT], self.Vh[ADJOINT])
dm = dl.TestFunction(self.Vh[PARAMETER])
res_form = self.varf_handler(u, m, p)
out.zero()
dl.assemble( dl.derivative(res_form, m, dm), tensor=out)
[docs] def setLinearizationPoint(self,x, gauss_newton_approx):
""" Set the values of the state and parameter
for the incremental forward and adjoint solvers. """
x_fun = [vector2Function(x[i], self.Vh[i]) for i in range(3)]
f_form = self.varf_handler(*x_fun)
g_form = [None,None,None]
for i in range(3):
g_form[i] = dl.derivative(f_form, x_fun[i])
self.A, dummy = dl.assemble_system(dl.derivative(g_form[ADJOINT],x_fun[STATE]), g_form[ADJOINT], self.bc0)
self.At, dummy = dl.assemble_system(dl.derivative(g_form[STATE],x_fun[ADJOINT]), g_form[STATE], self.bc0)
self.C = dl.assemble(dl.derivative(g_form[ADJOINT],x_fun[PARAMETER]))
[bc.zero(self.C) for bc in self.bc0]
if self.solver_fwd_inc is None:
self.solver_fwd_inc = self._createLUSolver()
self.solver_adj_inc = self._createLUSolver()
self.solver_fwd_inc.set_operator(self.A)
self.solver_adj_inc.set_operator(self.At)
if gauss_newton_approx:
self.Wuu = None
self.Wmu = None
self.Wmm = None
else:
self.Wuu = dl.assemble(dl.derivative(g_form[STATE],x_fun[STATE]))
[bc.zero(self.Wuu) for bc in self.bc0]
Wuu_t = Transpose(self.Wuu)
[bc.zero(Wuu_t) for bc in self.bc0]
self.Wuu = Transpose(Wuu_t)
self.Wmu = dl.assemble(dl.derivative(g_form[PARAMETER],x_fun[STATE]))
Wmu_t = Transpose(self.Wmu)
[bc.zero(Wmu_t) for bc in self.bc0]
self.Wmu = Transpose(Wmu_t)
self.Wmm = dl.assemble(dl.derivative(g_form[PARAMETER],x_fun[PARAMETER]))
[docs] def solveIncremental(self, out, rhs, is_adj):
""" If :code:`is_adj == False`:
Solve the forward incremental system:
Given :math:`u, m`, find :math:`\\tilde{u}` such that
.. math:: \\delta_{pu} F(u, m, p; \\hat{p}, \\tilde{u}) = \\mbox{rhs},\\quad \\forall \\hat{p}.
If :code:`is_adj == True`:
Solve the adjoint incremental system:
Given :math:`u, m`, find :math:`\\tilde{p}` such that
.. math:: \\delta_{up} F(u, m, p; \\hat{u}, \\tilde{p}) = \\mbox{rhs},\\quad \\forall \\hat{u}.
"""
if is_adj:
self.n_calls["incremental_adjoint"] += 1
self.solver_adj_inc.solve(out, rhs)
else:
self.n_calls["incremental_forward"] += 1
self.solver_fwd_inc.solve(out, rhs)
[docs] def apply_ij(self,i,j, dir, out):
"""
Given :math:`u, m, p`; compute
:math:`\\delta_{ij} F(u, m, p; \\hat{i}, \\tilde{j})` in the direction :math:`\\tilde{j} =` :code:`dir`,
:math:`\\forall \\hat{i}`.
"""
KKT = {}
KKT[STATE,STATE] = self.Wuu
KKT[PARAMETER, STATE] = self.Wmu
KKT[PARAMETER, PARAMETER] = self.Wmm
KKT[ADJOINT, STATE] = self.A
KKT[ADJOINT, PARAMETER] = self.C
if i >= j:
if KKT[i,j] is None:
out.zero()
else:
KKT[i,j].mult(dir, out)
else:
if KKT[j,i] is None:
out.zero()
else:
KKT[j,i].transpmult(dir, out)
[docs] def apply_ijk(self,i,j,k, x, jdir, kdir, out):
x_fun = [vector2Function(x[ii], self.Vh[ii]) for ii in range(3)]
idir_fun = dl.TestFunction(self.Vh[i])
jdir_fun = vector2Function(jdir, self.Vh[j])
kdir_fun = vector2Function(kdir, self.Vh[k])
res_form = self.varf_handler(*x_fun)
form = dl.derivative(
dl.derivative(
dl.derivative(res_form, x_fun[i], idir_fun),
x_fun[j], jdir_fun),
x_fun[k], kdir_fun)
out.zero()
dl.assemble(form, tensor=out)
if i in [STATE,ADJOINT]:
[bc.apply(out) for bc in self.bc0]
[docs] def _createLUSolver(self):
return PETScLUSolver(self.Vh[STATE].mesh().mpi_comm() )