Minor Release

examples/meshes/create_structured_hohlraum_2d_mesh.py

@@ -0,0 +1,137 @@
+import pygmsh as pg
+import numpy as np
+import itertools
+import os
+from optparse import OptionParser
+
+
+def add_block(x0, y0, x_len, y_len, char_length, geom):
+    coords = np.array([
+        [x0, y0, 0.0],
+        [x0 + x_len, y0, 0.0],
+        [x0 + x_len, y0 + y_len, 0.0],
+        [x0, y0 + y_len, 0.0]
+    ])
+    return geom.add_polygon(coords, char_length)
+
+
+def main():
+    print("---------- Start Creating the mesh ------------")
+    print("Parsing options")
+    # --- parse options ---
+    parser = OptionParser()
+    parser.add_option("-o", "--output_name", dest="output_name", default="hohlraum")
+    parser.add_option("-c", "--char_length", dest="char_length", default=0.025)
+    (options, args) = parser.parse_args()
+
+    options.output_name = str(options.output_name)
+    options.char_length = float(options.char_length)
+    char_length = options.char_length
+    geom = pg.opencascade.Geometry()
+    domain = add_block(-0.1, -0.05, 1.45, 1.4, char_length, geom)
+
+    boxes = [domain]
+    lines = []
+
+    # add interior polygon
+    coords = np.array([
+        [0.85, 0.25, 0.0],
+        [0.45, 0.25, 0.0],
+        [0.45, 1.05, 0.0],
+        [0.85, 1.05, 0.0],
+        [0.85, 1.0, 0.0],
+        [0.5, 1.0, 0.0],
+        [0.5, 0.3, 0.0],
+        [0.85, 0.3, 0.0],
+    ])
+    boxes.append(geom.add_polygon(coords, char_length))
+
+    # add outer polygon
+    coords = np.array([
+        [0.0, 0.0, 0.0],
+        [1.3, 0.0, 0.0],
+        [1.3, 1.3, 0.0],
+        [0.0, 1.3, 0.0],
+        [0.0, 1.25, 0.0],
+        [1.25, 1.25, 0.0],
+        [1.25, 0.05, 0.0],
+        [0.0, 0.05, 0.0],
+    ])
+    boxes.append(geom.add_polygon(coords, char_length))
+
+    """
+    points = []
+    points.append(geom.add_point([1.0, 1.05, 0.0], lcar=char_length))
+    points.append(geom.add_point([1.05, 1.05, 0.0], lcar=char_length))
+    points.append(geom.add_point([1.05, 1.0, 0.0], lcar=char_length))
+    points.append(geom.add_point([1.0, 1.0, 0.0], lcar=char_length))
+    lines = []
+    lines.append(geom.add_line(points[0], points[1]))
+    lines.append(geom.add_line(points[1], points[2]))
+    lines.append(geom.add_line(points[2], points[3]))
+    lines.append(geom.add_line(points[3], points[0]))
+
+    lineloop = geom.add_line_loop(lines)
+    surf = geom.add_surface(lineloop)
+    geom.add_physical(surf)
+    """
+    # points = []
+    # points.append(geom.add_point([0.0, 0.0, 0.0], lcar=char_length))
+    # points.append(geom.add_point([1.3, 0.0, 0.0], lcar=char_length))
+    # points.append(geom.add_point([1.3, 1.3, 0.0], lcar=char_length))
+    # points.append(geom.add_point([0.0, 1.3, 0.0], lcar=char_length))
+    # points.append(geom.add_point([0.0, 1.25, 0.0], lcar=char_length))
+    # points.append(geom.add_point([0.05, 1.25, 0.0], lcar=char_length))
+    # points.append(geom.add_point([1.25, 1.25, 0.0], lcar=char_length))
+    # points.append(geom.add_point([1.25, 1.25, 0.0], lcar=char_length))
+    # points.append(geom.add_point([1.25, 0.05, 0.0], lcar=char_length))
+    # points.append(geom.add_point([0.05, 0.05, 0.0], lcar=char_length))
+    # points.append(geom.add_point([0.0, 0.05, 0.0], lcar=char_length))
+    #
+    # for i in range(len(points) - 1):
+    #    lines.append(geom.add_line(points[i], points[i + 1]))
+    # lines.append(geom.add_line(points[len(points) - 1], points[0]))
+    #
+    # lineloop = geom.add_line_loop(lines)
+    # surf = geom.add_surface(lineloop)
+    # geom.add_physical(surf)
+
+    # add left side absorption region
+    boxes.append(add_block(0, 0.25, 0.05, 0.8, char_length, geom))
+    boxes.append(add_block(-0.05, 0.05, 0.05, 0.2, char_length, geom))
+    boxes.append(add_block(-0.05, 1.05, 0.05, 0.2, char_length, geom))
+
+    # add line that denotes the inflow ghost cells
+    # coords = np.array([
+    #    [-0.01, -0.01, 0.0],
+    #    [-0.01, 1.31, 0.0]])
+
+    # point1 = geom.add_point((0.0, 0.05, 0.0), lcar=char_length)
+    # point2 = geom.add_point((0.0, 0.25, 0.0), lcar=char_length)
+    # t = geom.add_line(point1, point2)
+    # geom.add_physical_line([t], label="void")
+
+    # lines.append(t)
+
+    geom.boolean_fragments(boxes, [])
+    geom.add_physical(domain.lines, label="void")
+
+    # domain.lines.append(t)
+    # geom.add_raw_code("Transfinite Surface{:};")
+    # geom.add_raw_code('Mesh.Algorithm= 3;')
+    # geom.add_raw_code("Recombine Surface {:} = 0;")
+
+    geom.add_raw_code("Transfinite Surface{:};")
+    geom.add_raw_code('Mesh.Algorithm= 3;')
+    geom.add_raw_code("Recombine Surface {:} = 0;")
+
+    mesh_code = geom.get_code()
+    with open(options.output_name + ".geo", "w") as mesh_file:
+        mesh_file.write(mesh_code)
+    os.system('gmsh ' + options.output_name + '.geo -2 -format su2 -save_all')
+    # os.system('rm ' + options.output_name + '.geo')
+    print("---------- Successfully created the mesh ------------")
+
+
+if __name__ == '__main__':
+    main()

include/common/globalconstants.hpp

@@ -149,10 +149,14 @@ inline std::map<std::string, ENTROPY_NAME> Entropy_Map{
     { "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLTZMANN", MAXWELL_BOLTZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
 
 // Optimizer
-enum OPTIMIZER_NAME { NEWTON, REGULARIZED_NEWTON, PART_REGULARIZED_NEWTON, ML };
-
-inline std::map<std::string, OPTIMIZER_NAME> Optimizer_Map{
-    { "NEWTON", NEWTON }, { "REGULARIZED_NEWTON", REGULARIZED_NEWTON }, { "PART_REGULARIZED_NEWTON", PART_REGULARIZED_NEWTON }, { "ML", ML } };
+enum OPTIMIZER_NAME { NEWTON, REGULARIZED_NEWTON, PART_REGULARIZED_NEWTON, REDUCED_NEWTON, REDUCED_PART_REGULARIZED_NEWTON, ML };
+
+inline std::map<std::string, OPTIMIZER_NAME> Optimizer_Map{ { "NEWTON", NEWTON },
+                                                            { "REGULARIZED_NEWTON", REGULARIZED_NEWTON },
+                                                            { "PART_REGULARIZED_NEWTON", PART_REGULARIZED_NEWTON },
+                                                            { "REDUCED_NEWTON", REDUCED_NEWTON },
+                                                            { "REDUCED_PART_REGULARIZED_NEWTON", REDUCED_PART_REGULARIZED_NEWTON },
+                                                            { "ML", ML } };
 
 // Volume output
 enum VOLUME_OUTPUT { ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS, MEDICAL };

include/datagenerator/datageneratorbase.hpp

@@ -58,6 +58,8 @@ class DataGeneratorBase
     NewtonOptimizer* _optimizer; /*!< @brief Class to solve minimal entropy problem */
     EntropyBase* _entropy;       /*!< @brief Class to handle entropy functional evaluations */
 
+    bool _reducedSampling; /*!< @brief Flag to show if the reduced optimizer is used */
+
     // Main methods
     virtual void SampleMultiplierAlpha(); /*!< @brief Sample Lagrange multipliers alpha */
     void ComputeRealizableSolution();     /*!< @brief make u the realizable moment to alpha, since Newton has roundoff errors. */
@@ -70,5 +72,7 @@ class DataGeneratorBase
     virtual void ComputeMoments() = 0;           /*!< @brief Pre-Compute Moments at all quadrature points. */
     bool ComputeEVRejection( unsigned idx_set ); /*!< @brief Evalute rejection criterion based on  the smallest Eigenvalue of the Hessian
                                                     corresponding to alpha[idx_set]. */
+    bool ComputeReducedEVRejection( VectorVector& redMomentBasis, Vector& redAlpha ); /*!< @brief Evalute rejection criterion based on  the smallest
+                                                     Eigenvalue of the reduced Hessian corresponding to alpha[idx_set]. */
 };
 #endif    // DATAGENERATOR_H

include/optimizers/neuralnetworkoptimizer.hpp

@@ -17,7 +17,7 @@ class NeuralNetworkOptimizer : public OptimizerBase
 
     void Solve( Vector& alpha, Vector& u, const VectorVector& moments, unsigned idx_cell = 0 ) override;
 
-    void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments ) override;
+    void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments, Vector& alpha_norms ) override;
 
     /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
      *  @param sol moment vector
@@ -59,7 +59,7 @@ class NeuralNetworkOptimizer : public OptimizerBase
 
     inline void Solve( Vector& alpha, Vector& u, const VectorVector& moments, unsigned idx_cell = 0 ) override{};
 
-    inline void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments ) override{};
+    inline void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments, Vector& alpha_norms ) override{};
 
     /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
      *  @param sol moment vector

include/optimizers/newtonoptimizer.hpp

@@ -19,7 +19,7 @@ class NewtonOptimizer : public OptimizerBase
     ~NewtonOptimizer();
 
     void Solve( Vector& alpha, Vector& sol, const VectorVector& moments, unsigned idx_cell = 0 ) override;
-    void SolveMultiCell( VectorVector& alpha, VectorVector& sol, const VectorVector& moments ) override;
+    void SolveMultiCell( VectorVector& alpha, VectorVector& sol, const VectorVector& moments, Vector& alpha_norms ) override;
 
     /*! @brief Computes the objective function
                 grad = <eta(alpha*m)> - alpha*sol */

include/optimizers/optimizerbase.hpp

@@ -34,7 +34,7 @@ class OptimizerBase
      *  @param   alpha vector where the solution Lagrange multipliers are saved to.
      *  @param   u  moment vector
      *  @param   moments  VectorVector to the moment basis evaluated at all quadpoints    */
-    virtual void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments ) = 0;
+    virtual void SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& moments, Vector& alpha_norms ) = 0;
 
     /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
      *  @param sol moment vector

include/optimizers/partregularizednewtonoptimizer.hpp

@@ -30,7 +30,7 @@ class PartRegularizedNewtonOptimizer : public NewtonOptimizer
     /*! @brief In 1D, this function scales the quadrature weigths to compute the entropy integrals in arbitrary (bounded) intervals
         @param leftBound : left boundary of the interval
         @param rightBound : right boundary of the interval*/
-    void ScaleQuadWeights( double leftBound, double rightBound );
+    // void ScaleQuadWeights( double leftBound, double rightBound );
 
     /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
      *  @param sol moment vector

include/optimizers/reducednewtonoptimizer.hpp

@@ -0,0 +1,40 @@
+/*!
+ * @file newtonoptimizer.h
+ * @brief class for solving the minimal entropy optimization problem using a newton optimizer with line search.
+ * @author S. Schotthöfer
+ */
+
+#ifndef REDUCEDNEWTONOPTIMIZER_H
+#define REDUCEDNEWTONOPTIMIZER_H
+
+#include "newtonoptimizer.hpp"
+
+class ReducedNewtonOptimizer : public NewtonOptimizer
+{
+  public:
+    ReducedNewtonOptimizer( Config* settings );
+
+    ~ReducedNewtonOptimizer();
+
+    /*! @brief Computes the objective function
+                grad = <eta(alpha*m)> - alpha*sol */
+    virtual double ComputeObjFunc( Vector& alpha, Vector& sol, const VectorVector& moments ) override;
+
+    /*! @brief Computes hessian of objective function and stores it in hessian
+        grad = <mXm*eta*'(alpha*m)> */
+    virtual void ComputeHessian( Vector& alpha, const VectorVector& moments, Matrix& hessian ) override;
+
+    /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
+     *  @param sol moment vector
+     *  @param alpha Lagrange multipliers
+     *  @param moments Moment basis
+     */
+    virtual void ReconstructMoments( Vector& sol, const Vector& alpha, const VectorVector& moments ) override;
+
+  protected:
+    /*! @brief Computes gradient of objective function and stores it in grad
+                grad = <m*eta*'(alpha*m)> - sol */
+    virtual void ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) override;
+};
+
+#endif    // REDUCEDNEWTONOPTIMIZER_H

include/optimizers/reducedpartregularizednewtonoptimizer.hpp

@@ -0,0 +1,48 @@
+/*!
+ * @file partregularizednewtonoptimizer.h
+ * @brief class for solving the minimal entropy optimization problem using a partly regularized (order zero moment is not regularized) newton
+ * optimizer with line search.
+ * @author S. Schotthöfer
+ */
+
+#ifndef REDUCEDPARTREGULARIZEDNEWTONOPTIMIZER_H
+#define REDUCEDPARTREGULARIZEDNEWTONOPTIMIZER_H
+
+#include "reducednewtonoptimizer.hpp"
+
+class ReducedPartRegularizedNewtonOptimizer : public ReducedNewtonOptimizer
+{
+  public:
+    ReducedPartRegularizedNewtonOptimizer( Config* settings );
+
+    ~ReducedPartRegularizedNewtonOptimizer();
+
+    /*! @brief Computes the objective function
+                grad = <eta(alpha*m)> - alpha*sol  + gamma/2*norm(alpha)*/
+    double ComputeObjFunc( Vector& alpha, Vector& sol, const VectorVector& moments ) override;
+
+    /*! @brief Computes hessian of objective function and stores it in hessian
+        grad = <mXm*eta*'(alpha*m)> */
+    void ComputeHessian( Vector& alpha, const VectorVector& moments, Matrix& hessian ) override;
+
+    /*! @brief In 1D, this function scales the quadrature weigths to compute the entropy integrals in arbitrary (bounded) intervals
+        @param leftBound : left boundary of the interval
+        @param rightBound : right boundary of the interval*/
+    // void ScaleQuadWeights( double leftBound, double rightBound );
+
+    /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha
+     *  @param sol moment vector
+     *  @param alpha Lagrange multipliers
+     *  @param moments Moment basis
+     */
+    void ReconstructMoments( Vector& sol, const Vector& alpha, const VectorVector& moments ) override;
+
+  private:
+    /*! @brief Computes gradient of objective function and stores it in grad
+                grad = <m*eta*'(alpha*m)> - sol + _gamma*alpha */
+    void ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) override;
+
+    double _gamma; /*!<  @brief Regularization Parameter*/
+};
+
+#endif    // REDUCEDPARTREGULARIZEDNEWTONOPTIMIZER_H

src/common/config.cpp

@@ -301,8 +301,7 @@ void Config::SetConfigOptions() {
     /*! @brief Neural Model Gamma \n DESCRIPTION:  Specifies regularization parameter for neural networks \n DEFAULT 0 (values: 0,1,2,3)
     \ingroup Config */
     AddUnsignedShortOption( "NEURAL_MODEL_GAMMA", _neuralGamma, 0 );
-    /*! @brief Neural Model Gamma \n DESCRIPTION:  Specifies regularization parameter for neural networks \n DEFAULT 0 (values: 0,1,2,3)
-    \ingroup Config */
+    /*! @brief NEURAL_MODEL_ENFORCE_ROTATION_SYMMETRY \n DESCRIPTION:  Flag to enforce rotational symmetry \n DEFAULT false \ingroup Config */
     AddBoolOption( "NEURAL_MODEL_ENFORCE_ROTATION_SYMMETRY", _enforceNeuralRotationalSymmetry, false );
 
     // Mesh related options
@@ -543,6 +542,15 @@ void Config::SetPostprocessing() {
         delete quad;
     }
 
+    // Optimizer Postprocessing
+    {
+        if( ( _entropyOptimizerName == REDUCED_NEWTON || _entropyOptimizerName == REDUCED_PART_REGULARIZED_NEWTON ) &&
+            _entropyName != MAXWELL_BOLTZMANN ) {
+            std::string msg = "Reduction of the optimization problen only possible with Maxwell Boltzmann Entropy" + std::to_string( _dim ) + ".";
+            ErrorMessages::Error( msg, CURRENT_FUNCTION );
+        }
+    }
+
     // --- Solver setup ---
     {
         if( GetSolverName() == PN_SOLVER && GetSphericalBasisName() != SPHERICAL_HARMONICS ) {
@@ -784,6 +792,19 @@ void Config::SetPostprocessing() {
         if( _regularizerGamma <= 0.0 ) {
             ErrorMessages::Error( "REGULARIZER_GAMMA must be positive.", CURRENT_FUNCTION );
         }
+
+        if( _entropyOptimizerName == ML ) {
+            // set regularizer gamma to the correct value
+            switch( _neuralGamma ) {
+                case 0: _regularizerGamma = 0; break;
+                case 1: _regularizerGamma = 0.1; break;
+                case 2: _regularizerGamma = 0.01; break;
+                case 3: _regularizerGamma = 0.001; break;
+            }
+        }
+        if( _entropyOptimizerName == NEWTON ) {
+            _regularizerGamma = 0.0;
+        }
     }
 }
 

src/common/mesh.cpp

@@ -17,7 +17,6 @@ Mesh::Mesh( std::vector<Vector> nodes,
     ComputeCellMidpoints();
     ComputeConnectivity();
     ComputeBounds();
-    // ComputeCellInterfaceMidpoints();
 }
 
 Mesh::~Mesh() {}
@@ -58,7 +57,7 @@ void Mesh::ComputeConnectivity() {
     }
 
     // determine neighbor cells and normals with MPI and OpenMP
-    #pragma omp parallel for
+#pragma omp parallel for
     for( unsigned i = mpiCellStart; i < mpiCellEnd; ++i ) {
         std::vector<unsigned>* cellsI = &sortedCells[i];
         unsigned ctr                  = 0;