Renaming functions and classes

.github/workflows/run_examples.sh

@@ -5,8 +5,7 @@ export HIPPYLIB_PATH=$(pwd)/hippylib
 # Run poisson example
 cd examples/poisson
 
-python3 driver_poisson_mean_serial.py
+python3 driver_poisson_mean.py
 
-mpirun -n 2 python3 driver_poisson_mean_mpi.py
+mpirun -n 2 python3 driver_poisson_mean.py
 
-python3 compare_results.py

.github/workflows/run_tests.sh

@@ -10,6 +10,6 @@ cd soupy/test
 $PYTHON -m unittest discover -v
 $PYTHON -m unittest discover -v -p 'ptest_*'
 
-mpirun -n 4 $PYTHON ptest_meanVarRiskMeasureSAA_MPI.py
+mpirun -n 4 $PYTHON ptest_meanVarRiskMeasureSAA.py
 mpirun -n 4 $PYTHON ptest_scipyCostWrapper.py
-mpirun -n 4 $PYTHON ptest_superquantileSAA_MPI.py
+mpirun -n 4 $PYTHON ptest_superquantileSAA.py

doc/soupy.modeling.rst

@@ -62,10 +62,10 @@ soupy.modeling.riskMeasure
     :show-inheritance:
 
 
-soupy.modeling.meanVarRiskMeasure
----------------------------------------
+soupy.modeling.meanVarRiskMeasureStochastic
+--------------------------------------------
 
-.. automodule:: soupy.modeling.meanVarRiskMeasure
+.. automodule:: soupy.modeling.meanVarRiskMeasureStochastic
     :members:
     :undoc-members:
     :show-inheritance:

examples/hyperelasticity/driver_hyperelasticity_deterministic.py → examples/hyperelasticity/driver_hyperelasticity.py

@@ -9,20 +9,43 @@
 #
 # SOUPy 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 3.0 dated June 1991.
+# Software Foundation) version 3.0 dated June 2007.
 
 
 """
-This example implements a deterministic optimization problem for material 
-design in hyperelasticity. The objective is to minimize the compliance of
-the structure subject. The design variable takes values in [0, 1], which 
-selects between two different material properties. 
+This example implements an optimization problem for material 
+design of a hyperelastic beam. The objective is to minimize the compliance of
+the structure subject to uncertain loading. 
+
+The design variable takes values in [0, 1], which selects between 
+two different material properties. A nominal external load 
+is prescribed to be at the center of the beam, but is scaled by 
+a multiplicative Gaussian random field.
+
+- Set risk measure flag :code:`-r` to :code:`deterministic` 
+    for Deterministic optimization using the mean parameter value 
+
+- Set risk measure flag :code:`-r` to :code:`mean_var` for optimization of 
+    the mean + variance risk measure 
+
+- Flag :code:`beta` controls the variance weighting 
+
+- Flag :code:`penalization` controls the l2 penalization on the design variable
 
 See :code:`setupHyperelasticityProblem.py` for problem settings, including
 mesh, geometry, and solver properties. 
 
 The example uses a custom forward solver. See 
 :code:`hyperelasticityControlPDE.py` for the PDE definition and solver
+
+This example also shows how to use the :code:`ScipyCostWrapper` to convert
+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
+
 """
 
 
@@ -53,64 +76,91 @@
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
-    parser.add_argument('-n', '--maxiter', type=int, default=100, help="Maximum number of SD iterations")
+    parser.add_argument('-r', '--risk_measure', type=str, default="deterministic", choices=["deterministic", "mean_var"], help="Risk measure type")
+    parser.add_argument('-n', '--sample_size', type=int, default=32, help="Sample size for risk measure computation")
+    parser.add_argument('-b', '--beta', type=float, default=1.0, help="Variance weight for risk measure")
     parser.add_argument('-q', '--qoi_type', type=str, default="stiffness", choices=["all", "stiffness", "point"])
     parser.add_argument('-p', '--penalization', type=float, default=1e-1, help="Scaling of penalization")
+
+    parser.add_argument('--maxiter', type=int, default=100, help="Maximum number of SD iterations")
     parser.add_argument('--show', default=False, action="store_true", help="Show figure")
     args = parser.parse_args()
 
+    save_dir = "results_%s" %(args.risk_measure)
+    os.makedirs(save_dir, exist_ok=True)
+
     # Create mesh and setup model components 
-    comm_mesh = MPI.COMM_WORLD
+    comm_mesh = MPI.COMM_SELF
+    comm_sampler = MPI.COMM_WORLD 
 
     settings = hyperelasticity_problem_settings()
     settings["qoi_type"] = args.qoi_type
     mesh, Vh, hyperelasticity_varf, control_model, prior = setup_hyperelasticity_problem(settings, comm_mesh)
 
     l2_penalty = soupy.L2Penalization(Vh, args.penalization)
-    pde_cost = soupy.DeterministicControlCostFunctional(control_model, prior, l2_penalty)
 
+    if args.risk_measure == "deterministic":
+        pde_cost = soupy.DeterministicControlCostFunctional(control_model, prior, l2_penalty)
+
+    else:
+        # Use the mean variance risk measure to assemble cost 
+        risk_settings = soupy.meanVarRiskMeasureSAASettings()
+        risk_settings["beta"] = args.beta
+        risk_settings["sample_size"] = args.sample_size
+        risk_measure = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings, comm_sampler=comm_sampler)
+        pde_cost = soupy.RiskMeasureControlCostFunctional(risk_measure, l2_penalty)
+
+    # ------------------  Using the scipy cost interface ------------------ # 
     scipy_cost = soupy.ScipyCostWrapper(pde_cost, verbose=True)
     box_bounds = scipy.optimize.Bounds(lb=0.0, ub=1.0)
 
+    # Solve the PDE for the beam using entirely soft material 
     x = control_model.generate_vector()
     x[soupy.PARAMETER].axpy(1.0, prior.mean)
     control_model.solveFwd(x[soupy.STATE], x)
     disp_fun_init = hp.vector2Function(x[soupy.STATE], Vh[soupy.STATE])
 
+    # Solve optimal design problem using 0.5 as initial guess 
     z0_np = x[soupy.CONTROL].get_local() + 0.5
     options = {'gtol' : 1e-12, 'disp' : True, 'maxiter' : args.maxiter}
     results = scipy.optimize.minimize(scipy_cost.function(), z0_np, 
             method="L-BFGS-B", jac=scipy_cost.jac(), 
             bounds=box_bounds, options=options)
+
+    # Optimal design 
     z_opt = results['x']     
 
+    # --------------------------------------------------------------------- # 
+
     x[soupy.CONTROL].set_local(z_opt)
     x[soupy.CONTROL].apply("")
 
     control_model.solveFwd(x[soupy.STATE], x)
     z_fun = dl.Function(Vh[soupy.CONTROL], x[soupy.CONTROL])
-
-    with dl.HDF5File(mesh.mpi_comm(), "z_opt.h5", "w") as save_file:
-        save_file.write(z_fun, "control")
+
 
     # ---------- postprocessing the plotting ----------- # 
+
     plt.figure()
     ax = dl.plot(disp_fun_init, mode="displacement")
     plt.colorbar(ax)
     plt.title("Displacement with soft material")
-    plt.savefig("u_soft.png")
+    plt.savefig("%s/u_soft.png" %(save_dir))
 
     plt.figure()
     ax = dl.plot(hp.vector2Function(x[soupy.CONTROL], Vh[soupy.CONTROL]))
     plt.colorbar(ax)
     plt.title("Optimal design")
-    plt.savefig("z_opt.png")
+    plt.savefig("%s/z_opt.png" %(save_dir))
 
     plt.figure()
     ax = dl.plot(hp.vector2Function(x[soupy.STATE], Vh[soupy.STATE]), mode="displacement")
     plt.colorbar(ax)
     plt.title("Optimal displacement")
-    plt.savefig("u_opt.png")
+    plt.savefig("%s/u_opt.png" %(save_dir))
+
+    if comm_sampler.Get_rank() == 0:
+        np.save('%s/z_opt.npy' %(save_dir), z_opt)
 
     if args.show:
         plt.show()

examples/hyperelasticity/hyperelasticityControlPDE.py

@@ -9,7 +9,7 @@
 #
 # SOUPy 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 3.0 dated June 1991.
+# Software Foundation) version 3.0 dated June 2007.
 
 import pickle
 import sys, os

examples/hyperelasticity/setupHyperelasticityProblem.py

@@ -9,7 +9,7 @@
 #
 # SOUPy 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 3.0 dated June 1991.
+# Software Foundation) version 3.0 dated June 2007.
 
 import pickle
 import sys, os
@@ -67,7 +67,7 @@ def hyperelasticity_problem_settings():
     settings["geometry"] = {"lx" : 2.0, "ly" : 0.5, "lz" : 0.25, "dim" : 2}
     settings["mesh"] = {"nx" : 96, "ny" : 24, "nz" : 12}
     settings["solver"] = {"backtrack" : True, "load_steps" : [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 0.9375, 1.0]}
-    settings["uncertainty"] = {"gamma" : 0.1, "delta" : 1, "robin_bc" : True}
+    settings["uncertainty"] = {"gamma" : 0.2, "delta" : 1.0, "robin_bc" : True}
     return settings
 
 

examples/poisson/driver_poisson_mean_mpi.py → examples/poisson/driver_poisson_mean.py

@@ -9,14 +9,16 @@
 #
 # SOUPy 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 3.0 dated June 1991.
+# Software Foundation) version 3.0 dated June 2007.
 
 """
 Minimization of the mean risk measure using the SAA approximation for the 
 linear poisson problem with a log-normal conductivity field. 
 
-This driver uses MPI to parallelize the sampling of the parameter field.
-For a purely serial implementation, see :code:`driver_poisson_mean_serial.py`
+This driver supports MPI to parallelize the sampling of the parameter field.
+e.g. 
+
+mpirun -n 4 python driver_poisson_mean.py
 """
 
 import sys 
@@ -45,7 +47,7 @@
 SAMPLE_SIZE = 4
 PENALTY_WEIGHT = 1e-3
 
-RESULTS_DIRECTORY = "results_mpi"
+RESULTS_DIRECTORY = "results"
 
 # 0. MPI Communicators
 comm_mesh = MPI.COMM_SELF
@@ -93,7 +95,7 @@ def boundary(x, on_boundary):
 risk_settings = soupy.meanVarRiskMeasureSAASettings()
 risk_settings["beta"] = VARIANCE_WEIGHT
 risk_settings["sample_size"] = SAMPLE_SIZE 
-risk_measure = soupy.MeanVarRiskMeasureSAA_MPI(control_model, prior, risk_settings, comm_sampler=comm_sampler)
+risk_measure = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings, comm_sampler=comm_sampler)
 
 # 7. Define the penalization term 
 penalty = soupy.L2Penalization(Vh, PENALTY_WEIGHT)
@@ -103,7 +105,7 @@ def boundary(x, on_boundary):
 
 
 # 9. Define the optimizer 
-optimizer = soupy.BFGS(cost_functional)
+optimizer = soupy.InexactNewtonCG(cost_functional)
 
 # 10. Provide initial guess and solve 
 if comm_sampler.rank == 0 :