Updates to penalization and QoI classes

examples/hyperelasticity/driver_hyperelasticity.py

@@ -42,10 +42,16 @@
 a :code:`ControlCostFunctional` to be compatible with :code:`scipy.optimize`
 
 This driver supports MPI to parallelize the sampling of the parameter field.
-e.g. 
 
-mpirun -n 4 python driver_hyperelasticity.py
+To run deterministic: 
+python driver_hyperelasticity.py
 
+
+To run with mean + variance risk measure:
+python driver_hyperelasticity.py -r mean_var
+
+To run with mean + variance risk measure and parllel sampling (e.g.):
+mpirun -n 4 python driver_hyperelasticity.py -r mean_var
 """
 
 

examples/hyperelasticity/setupHyperelasticityProblem.py

@@ -89,13 +89,13 @@ def setup_qoi(mesh, Vh, settings):
     """
     if settings["qoi_type"] == "all":
         u0 = dl.Function(Vh[soupy.STATE]).vector()
-        qoi = soupy.L2MisfitControlQoI(mesh, Vh, u0)
+        qoi = soupy.L2MisfitControlQoI(Vh, u0)
     elif settings["qoi_type"] == "stiffness":
         stiffness = Stiffnessvarf()
-        qoi = soupy.VariationalControlQoI(mesh, Vh, stiffness)
+        qoi = soupy.VariationalControlQoI(Vh, stiffness)
     elif settings["qoi_type"] == "point":
         local_displacement = LocalDisplacementVarf(settings["obs"]["location"], settings["obs"]["width"])
-        qoi = soupy.VariationalControlQoI(mesh, Vh, local_displacement)
+        qoi = soupy.VariationalControlQoI(Vh, local_displacement)
     else:
         raise ValueError("Settings qoi type not available")
 

examples/poisson/driver_poisson_mean.py

@@ -86,7 +86,7 @@ def boundary(x, on_boundary):
 u_target = dl.Expression("x[1] + sin(k*x[0]) * sin(k*x[1])", k=1.5*np.pi, degree=2, mpi_comm=comm_mesh)
 u_target_function = dl.interpolate(u_target, Vh_STATE)
 u_target_vector = u_target_function.vector()
-qoi = soupy.L2MisfitControlQoI(mesh, Vh, u_target_vector)
+qoi = soupy.L2MisfitControlQoI(Vh, u_target_vector)
 
 # 5. Define the ControlModel
 control_model = soupy.ControlModel(pde, qoi)

examples/semilinear_elliptic/driver_semilinear_cvar.py

@@ -135,7 +135,7 @@ def print_on_root(print_str, mpi_comm=MPI.COMM_WORLD):
     print_on_root("Making target, QoI, and problem", comm_sampler)
     u_target_expr = get_target(args.target, args.param)
     u_target = dl.interpolate(u_target_expr, Vh[hp.STATE])
-    qoi = soupy.L2MisfitControlQoI(mesh, Vh, u_target.vector())
+    qoi = soupy.L2MisfitControlQoI(Vh, u_target.vector())
     control_model = soupy.ControlModel(pde, qoi)
 
     print_on_root("Making SAA risk measure and cost", comm_sampler)

soupy/modeling/__init__.py

@@ -29,7 +29,7 @@
 
 from .meanVarRiskMeasureSAA import meanVarRiskMeasureSAASettings, MeanVarRiskMeasureSAA 
 
-from .penalization import Penalization, L2Penalization, WeightedL2Penalization, MultiPenalization
+from .penalization import Penalization, L2Penalization, WeightedL2Penalization, MultiPenalization, VariationalPenalization
 
 from .riskMeasure import RiskMeasure
 

soupy/modeling/controlQoI.py

@@ -100,19 +100,17 @@ class VariationalControlQoI(ControlQoI):
     """
     Class for a QoI defined by its variational form 
     """
-    def __init__(self, mesh, Vh, form_handler):
+    def __init__(self, Vh, form_handler):
         """
         Constructor
 
-        :param mesh: The mesh object
         :param Vh: List of function spaces for the state, parameter, 
             adjoint, and optimization variables
         :type Vh: list of :py:class:`dolfin.FunctionSpace`
         :param form_handler: The form handler for the variational form with a 
             :code:`__call__` method that takes as input the state, parameter, and control variables
             as functions and returns the variational form 
         """
-        self.mesh = mesh
         self.Vh = Vh
         self.x = [dl.Function(Vh[STATE]).vector(), dl.Function(Vh[PARAMETER]).vector(),
                   dl.Function(Vh[ADJOINT]).vector(), dl.Function(Vh[CONTROL]).vector()]
@@ -254,18 +252,16 @@ class L2MisfitControlQoI(ControlQoI):
 
     where :math:`u_d` is the reference state.
     """
-    def __init__(self, mesh, Vh, ud):
+    def __init__(self, Vh, ud):
         """
         Constructor
 
-        :param mesh: The mesh object
         :param Vh: List of function spaces for the state, parameter, 
             adjoint, and optimization variables
         :type Vh: list of :py:class:`dolfin.FunctionSpace`
         :param ud: The reference state as a vector
         :type ud: :py:class:`dolfin.Vector`
         """
-        self.mesh = mesh
         self.Vh = Vh
         self.x = [dl.Function(Vh[STATE]).vector(), dl.Function(Vh[PARAMETER]).vector(),
                   dl.Function(Vh[ADJOINT]).vector(), dl.Function(Vh[CONTROL]).vector()]

soupy/modeling/penalization.py

@@ -33,6 +33,80 @@ def hessian(self, z, zhat, out):
         raise NotImplementedError("Child class should implement hessian")
 
 
+
+class VariationalPenalization(Penalization):
+    """
+    Penalization given by a variational form in terms of the control variable :math:`z`
+    """
+    def __init__(self, Vh, form_handler):
+        """
+        Constructor: 
+
+        :param Vh: List of function spaces for the state, parameter, 
+            adjoint, and optimization variables
+        :type Vh: list of :py:class:`dolfin.FunctionSpace`
+        :param form_handler: The form handler for the penalization with a
+            :code:`__call__` method that takes the control variable :code:`z` as a :py:class:`dolfin.Function` and returns the variational form. 
+        """
+
+        self.Vh = Vh 
+        self.form_handler = form_handler
+
+        self.z_fun = dl.Function(self.Vh[CONTROL])
+        self.z = self.z_fun.vector()
+
+        self.dz_fun = dl.Function(self.Vh[CONTROL])
+        self.dz = self.dz_fun.vector()
+
+        self.z_trial = dl.TrialFunction(self.Vh[CONTROL])
+        self.z_test = dl.TestFunction(self.Vh[CONTROL])
+
+
+        self.penalization_form = self.form_handler(self.z_fun)
+        self.grad_form = dl.derivative(self.penalization_form, self.z_fun, self.z_test)
+
+    def init_vector(self, z):
+        if isinstance(z, AugmentedVector):
+            pass 
+        else:
+            z.init(self.z_fun.vector().local_range())
+
+    def cost(self, z):
+        if isinstance(z, AugmentedVector):
+            z = z.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+        return dl.assemble(self.penalization_form)
+
+    def grad(self, z, out):
+        if isinstance(z, AugmentedVector):
+            out.set_scalar(0)
+            z = z.get_vector()
+            out = out.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+        dl.assemble(self.grad_form, tensor=out)
+
+    def hessian(self, z, zhat, out):
+        if isinstance(z, AugmentedVector):
+            out.set_scalar(0)
+            z = z.get_vector()
+            zhat = zhat.get_vector()
+            out = out.get_vector()
+
+        self.z.zero()
+        self.z.axpy(1.0, z)
+
+        self.dz.zero()
+        self.dz.axpy(1.0, zhat)
+
+        hessian_action = dl.derivative(self.grad_form, self.z_fun, self.dz_fun)
+        dl.assemble(hessian_action, tensor=out)
+
+
+
 class MultiPenalization(Penalization):
     """
     Penalization term for the sum of individual penalties
@@ -61,6 +135,8 @@ def __init__(self, Vh, penalization_list, alpha_list=None):
         else:
             self.alpha_list = [1.0] * len(self.penalization_list)
 
+    def init_vector(self, z):
+        self.penalization_list[0].init_vector(z)
 
     def cost(self, z):
         """

soupy/test/ptest_meanVarRiskMeasureSAA.py

@@ -75,7 +75,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         risk = MeanVarRiskMeasureSAA(model, prior, comm_sampler=self.comm_sampler)
@@ -96,7 +96,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 

soupy/test/ptest_scipyCostWrapper.py

@@ -93,7 +93,7 @@ def _sample_control(self, z):
         self.collective.bcast(z, root=0)
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 

soupy/test/ptest_superquantileSAA.py

@@ -103,7 +103,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 
@@ -279,7 +279,7 @@ def computeSuperquantileValue(self, sample_size, beta):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = superquantileRiskMeasureSAASettings()

soupy/test/test_SAACostFunctional.py

@@ -68,7 +68,7 @@ def _setup_pde_and_distributions(self, is_fwd_linear):
         return pde, prior, control_dist
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 

soupy/test/test_controlCostFunctional.py

@@ -76,7 +76,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(2.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
         z = model.generate_vector(CONTROL)
@@ -102,7 +102,7 @@ def finiteDifferenceCheck(self, qoi_varf, is_fwd_linear=True):
         settings['LINEAR'] = is_fwd_linear
 
         pde, prior, control_dist = setupPoissonPDEProblem(self.Vh, settings) 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
 
@@ -196,7 +196,7 @@ def testSavedSolution(self):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         cost = DeterministicControlCostFunctional(model, prior, penalization=None)
 

soupy/test/test_controlQoI.py

@@ -58,7 +58,7 @@ def finiteDifferenceGradient(self, nonlinear_qoi, i):
 
         # ud = dl.Function(self.Vh[STATE])
         # ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
-        # nonlinear_qoi = VariationalControlQoI(self.mesh, self.Vh, varf)
+        # nonlinear_qoi = VariationalControlQoI(self.Vh, varf)
 
 
         u_fun = dl.Function(self.Vh[STATE])
@@ -137,7 +137,7 @@ def testFiniteDifference(self):
         ud = dl.Function(self.Vh[STATE])
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
         for varf in varfs:
-            nonlinear_qoi = VariationalControlQoI(self.mesh, self.Vh, varf)
+            nonlinear_qoi = VariationalControlQoI(self.Vh, varf)
             for i in [STATE, PARAMETER, CONTROL]:
                 self.finiteDifferenceGradient(nonlinear_qoi, i)
                 for j in [STATE, PARAMETER, CONTROL]:
@@ -150,7 +150,7 @@ def testVariationalControlQoI(self):
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=5))
         chi = dl.Expression("0.5", degree=5)
         l2varf = L2MisfitVarfHandler(ud, chi=chi)
-        l2qoi = VariationalControlQoI(self.mesh, self.Vh, l2varf)
+        l2qoi = VariationalControlQoI(self.Vh, l2varf)
         u = dl.Function(self.Vh[STATE])
         m = dl.Function(self.Vh[PARAMETER])
         z = dl.Function(self.Vh[CONTROL])
@@ -166,7 +166,7 @@ def testVariationalControlQoI(self):
     def testL2MisfitQoI(self):
         ud = dl.Function(self.Vh[STATE])
         ud.interpolate(dl.Expression("sin(k*x[0])*cos(k*x[1])", k=math.pi, degree=2))
-        nonlinear_qoi = L2MisfitControlQoI(self.mesh, self.Vh, ud.vector())
+        nonlinear_qoi = L2MisfitControlQoI(self.Vh, ud.vector())
         for i in [STATE, PARAMETER, CONTROL]:
             self.finiteDifferenceGradient(nonlinear_qoi, i)
             for j in [STATE, PARAMETER, CONTROL]:

soupy/test/test_meanVarRiskMeasureStochastic.py

@@ -63,7 +63,7 @@ def testCostValue(self):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx 
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         risk = MeanVarRiskMeasureStochastic(model, prior)
@@ -91,7 +91,7 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = meanVarRiskMeasureStochasticSettings()

soupy/test/test_scipyCostWrapper.py

@@ -84,7 +84,7 @@ def _sample_control(self, z):
         self.control_prior.sample(self.control_noise, z)
 
     def _setup_control_model(self, pde, qoi_varf):
-        qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf)
+        qoi = VariationalControlQoI(self.Vh, qoi_varf)
         model = ControlModel(pde, qoi)
         return qoi, model
 

soupy/test/test_stochasticCostFunctional.py

@@ -65,7 +65,7 @@ def costValueCheck(self, sample_size, use_penalization):
         def l2norm(u,m,z):
             return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
         risk = MeanVarRiskMeasureStochastic(model, prior)
         alpha = 2.0
@@ -124,7 +124,7 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True):
         def l2norm(u,m,z):
             return u**2*dl.dx
 
-        qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm)
+        qoi = VariationalControlQoI(self.Vh, l2norm)
         model = ControlModel(pde, qoi)
 
         rm_settings = meanVarRiskMeasureStochasticSettings()