diff --git a/examples/meshes/create_structured_hohlraum_2d_mesh.py b/examples/meshes/create_structured_hohlraum_2d_mesh.py new file mode 100644 index 00000000..f4255445 --- /dev/null +++ b/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() diff --git a/include/common/globalconstants.hpp b/include/common/globalconstants.hpp index 8a76ede2..3304a669 100644 --- a/include/common/globalconstants.hpp +++ b/include/common/globalconstants.hpp @@ -149,10 +149,14 @@ inline std::map 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 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 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 }; diff --git a/include/datagenerator/datageneratorbase.hpp b/include/datagenerator/datageneratorbase.hpp index 7d0507bd..c069d9e6 100644 --- a/include/datagenerator/datageneratorbase.hpp +++ b/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 diff --git a/include/optimizers/neuralnetworkoptimizer.hpp b/include/optimizers/neuralnetworkoptimizer.hpp index 4a66390b..fd8cc99c 100644 --- a/include/optimizers/neuralnetworkoptimizer.hpp +++ b/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 diff --git a/include/optimizers/newtonoptimizer.hpp b/include/optimizers/newtonoptimizer.hpp index 400f7e37..5f15d40a 100644 --- a/include/optimizers/newtonoptimizer.hpp +++ b/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 = - alpha*sol */ diff --git a/include/optimizers/optimizerbase.hpp b/include/optimizers/optimizerbase.hpp index eca0fe76..54effb26 100644 --- a/include/optimizers/optimizerbase.hpp +++ b/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 diff --git a/include/optimizers/partregularizednewtonoptimizer.hpp b/include/optimizers/partregularizednewtonoptimizer.hpp index 027a80d6..5f7e8ea4 100644 --- a/include/optimizers/partregularizednewtonoptimizer.hpp +++ b/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 diff --git a/include/optimizers/reducednewtonoptimizer.hpp b/include/optimizers/reducednewtonoptimizer.hpp new file mode 100644 index 00000000..9926f157 --- /dev/null +++ b/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 = - 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 = */ + 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 = - sol */ + virtual void ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) override; +}; + +#endif // REDUCEDNEWTONOPTIMIZER_H diff --git a/include/optimizers/reducedpartregularizednewtonoptimizer.hpp b/include/optimizers/reducedpartregularizednewtonoptimizer.hpp new file mode 100644 index 00000000..0e9c6c09 --- /dev/null +++ b/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 = - 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 = */ + 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 = - sol + _gamma*alpha */ + void ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) override; + + double _gamma; /*!< @brief Regularization Parameter*/ +}; + +#endif // REDUCEDPARTREGULARIZEDNEWTONOPTIMIZER_H diff --git a/src/common/config.cpp b/src/common/config.cpp index 0ea69306..1ad40dd6 100644 --- a/src/common/config.cpp +++ b/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; + } } } diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index 3872d373..1e43d967 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -17,7 +17,6 @@ Mesh::Mesh( std::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* cellsI = &sortedCells[i]; unsigned ctr = 0; diff --git a/src/datagenerator/datageneratorbase.cpp b/src/datagenerator/datageneratorbase.cpp index c2a283ef..29f150ec 100644 --- a/src/datagenerator/datageneratorbase.cpp +++ b/src/datagenerator/datageneratorbase.cpp @@ -14,6 +14,8 @@ #include "entropies/entropybase.hpp" #include "optimizers/newtonoptimizer.hpp" #include "optimizers/partregularizednewtonoptimizer.hpp" +#include "optimizers/reducednewtonoptimizer.hpp" +#include "optimizers/reducedpartregularizednewtonoptimizer.hpp" #include "optimizers/regularizednewtonoptimizer.hpp" #include "quadratures/quadraturebase.hpp" #include "toolboxes/errormessages.hpp" @@ -62,10 +64,19 @@ DataGeneratorBase::DataGeneratorBase( Config* settings ) { _momentBasis = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) ); // Optimizer + _reducedSampling = false; switch( settings->GetOptimizerName() ) { case NEWTON: _optimizer = new NewtonOptimizer( settings ); break; case REGULARIZED_NEWTON: _optimizer = new RegularizedNewtonOptimizer( settings ); break; case PART_REGULARIZED_NEWTON: _optimizer = new PartRegularizedNewtonOptimizer( settings ); break; + case REDUCED_NEWTON: + _optimizer = new ReducedNewtonOptimizer( settings ); + _reducedSampling = true; + break; + case REDUCED_PART_REGULARIZED_NEWTON: + _optimizer = new ReducedPartRegularizedNewtonOptimizer( settings ); + _reducedSampling = true; + break; default: ErrorMessages::Error( "Optimizer choice not feasible for datagenerator.", CURRENT_FUNCTION ); break; } // Entropy @@ -126,6 +137,7 @@ void DataGeneratorBase::SampleMultiplierAlpha() { Vector alphaRed = Vector( _nTotalEntries - 1, 0.0 ); // local reduced alpha bool accepted = false; + bool normAccepted = false; while( !accepted ) { // Sample random multivariate uniformly distributed alpha between minAlpha and MaxAlpha. for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { @@ -137,6 +149,10 @@ void DataGeneratorBase::SampleMultiplierAlpha() { if( alphaRed[idx_sys - 1] < -1 * maxAlphaValue ) alphaRed[idx_sys - 1] = -1 * maxAlphaValue; } } + normAccepted = true; + if( _settings->GetUniformSamlping() ) { + if( norm( alphaRed ) > maxAlphaValue ) normAccepted = false; + } // Compute alpha_0 = log() // for maxwell boltzmann! only double integral = 0.0; // Integrate @@ -151,96 +167,23 @@ void DataGeneratorBase::SampleMultiplierAlpha() { } // Compute rejection criteria - accepted = ComputeEVRejection( idx_set ); - } - } - } - else { - // non normalized sampling - // TODO - ErrorMessages::Error( "Non-Normalized Alpha Sampling is not yet implemented.", CURRENT_FUNCTION ); - } - - // Cubeoid - /* - double minAlphaValue = -1 * maxAlphaValue; - // double maxAlphaValue = maxAlphaValue; - if( _settings->GetNormalizedSampling() ) { - // compute reduced version of alpha and m - if( _maxPolyDegree == 0 ) { - ErrorMessages::Error( "Normalized sampling not meaningful for M0 closure", CURRENT_FUNCTION ); - } - - VectorVector alphaRed = VectorVector( _setSize, Vector( _nTotalEntries - 1, 0.0 ) ); - VectorVector momentsRed = VectorVector( _nq, Vector( _nTotalEntries - 1, 0.0 ) ); - - for( unsigned idx_nq = 0; idx_nq < _nq; idx_nq++ ) { // copy (reduced) moments - for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { - momentsRed[idx_nq][idx_sys - 1] = _momentBasis[idx_nq][idx_sys]; - } - } - - // Sample alphaRed as uniform grid from [minAlphaValue, maxAlphaValue], then compute alpha_0 s.t. u_0 = 1 - double _gridSize = floor( pow( (double)_setSize, 1 / (double)( _nTotalEntries - 1 ) ) ); - double dalpha = ( maxAlphaValue - minAlphaValue ) / (double)_gridSize; - unsigned count = 0; - - switch( _nTotalEntries - 1 ) { - case 1: - for( unsigned idx_set = 0; idx_set < _gridSize; idx_set++ ) { - alphaRed[idx_set][0] = minAlphaValue + idx_set * dalpha; - } - break; - case 2: - count = 0; - for( unsigned i1 = 0; i1 < _gridSize; i1++ ) { - double alpha0 = minAlphaValue + i1 * dalpha; - for( unsigned i2 = 0; i2 < _gridSize; i2++ ) { - alphaRed[count][0] = alpha0; - alphaRed[count][1] = minAlphaValue + i2 * dalpha; - count++; + accepted = normAccepted; + if( normAccepted ) { + if( _reducedSampling ) { + accepted = ComputeReducedEVRejection( momentsRed, alphaRed ); } - } - break; - case 3: - count = 0; - for( unsigned i1 = 0; i1 < _gridSize; i1++ ) { - double alpha0 = minAlphaValue + i1 * dalpha; - for( unsigned i2 = 0; i2 < _gridSize; i2++ ) { - double alpha1 = minAlphaValue + i2 * dalpha; - for( unsigned i3 = 0; i3 < _gridSize; i3++ ) { - alphaRed[count][0] = alpha0; - alphaRed[count][1] = alpha1; - alphaRed[count][2] = minAlphaValue + i3 * dalpha; - count++; - } + else { + accepted = ComputeEVRejection( idx_set ); } } - break; - default: ErrorMessages::Error( "Not yet implemented!", CURRENT_FUNCTION ); - } - - // Compute alpha_0 = log() // for maxwell boltzmann! only - for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { - double integral = 0.0; - // Integrate (eta(eta'_*(alpha*m)) - for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { - integral += _entropy->EntropyPrimeDual( dot( alphaRed[idx_set], momentsRed[idx_quad] ) ) * _weights[idx_quad]; - } - _alpha[idx_set][0] = -log( integral ); // log trafo - - // copy all other alphas to the member - for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { - _alpha[idx_set][idx_sys] = alphaRed[idx_set][idx_sys - 1]; } } } else { // non normalized sampling // TODO - ErrorMessages::Error( "Not yet implemented!", CURRENT_FUNCTION ); + ErrorMessages::Error( "Non-Normalized Alpha Sampling is not yet implemented.", CURRENT_FUNCTION ); } - */ } void DataGeneratorBase::PrintLoadScreen() { @@ -250,28 +193,58 @@ void DataGeneratorBase::PrintLoadScreen() { } bool DataGeneratorBase::ComputeEVRejection( unsigned idx_set ) { - Matrix hessian = Matrix( _nTotalEntries, _nTotalEntries, 0.0 ); - // std::cout << idx_set << "\n"; - // std::cout << _alpha[idx_set] << "\n"; - // TextProcessingToolbox::PrintVectorVector( _momentBasis ); _optimizer->ComputeHessian( _alpha[idx_set], _momentBasis, hessian ); - // TextProcessingToolbox::PrintMatrix( hessian ); SymMatrix hessianSym( hessian ); // Bad solution, rewrite with less memory need Vector ew = Vector( _nTotalEntries, 0.0 ); eigen( hessianSym, ew ); if( min( ew ) < _settings->GetMinimalEVBound() ) { - // std::cout << "Sampling not accepted with EV:" << min( ew ) << std::endl; - // std::cout << hessianSym << std::endl; - // std::cout << "Current minimal accepted EV:" << _settings->GetMinimalEVBound() << std::endl; + return false; + } + return true; +} + +bool DataGeneratorBase::ComputeReducedEVRejection( VectorVector& redMomentBasis, Vector& redAlpha ) { + Matrix hessian = Matrix( _nTotalEntries - 1, _nTotalEntries - 1, 0.0 ); + _optimizer->ComputeHessian( redAlpha, redMomentBasis, hessian ); + SymMatrix hessianSym( hessian ); // Bad solution, rewrite with less memory need + Vector ew = Vector( _nTotalEntries - 1, 0.0 ); + eigen( hessianSym, ew ); + if( min( ew ) < _settings->GetMinimalEVBound() ) { return false; } return true; } void DataGeneratorBase::ComputeRealizableSolution() { + + if( _reducedSampling ) { + VectorVector momentsRed = VectorVector( _nq, Vector( _nTotalEntries - 1, 0.0 ) ); + + for( unsigned idx_nq = 0; idx_nq < _nq; idx_nq++ ) { // copy (reduced) moments + for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { + momentsRed[idx_nq][idx_sys - 1] = _momentBasis[idx_nq][idx_sys]; + } + } #pragma omp parallel for schedule( guided ) - for( unsigned idx_sol = 0; idx_sol < _setSize; idx_sol++ ) { - _optimizer->ReconstructMoments( _uSol[idx_sol], _alpha[idx_sol], _momentBasis ); + for( unsigned idx_sol = 0; idx_sol < _setSize; idx_sol++ ) { + Vector uSolRed( _nTotalEntries - 1, 0.0 ); + Vector alphaRed( _nTotalEntries - 1, 0.0 ); + for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { + alphaRed[idx_sys - 1] = _alpha[idx_sol][idx_sys]; + } + _optimizer->ReconstructMoments( uSolRed, alphaRed, momentsRed ); + + for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { + _uSol[idx_sol][idx_sys] = uSolRed[idx_sys - 1]; + } + _uSol[idx_sol][0] = 1.0; + } + } + else { +#pragma omp parallel for schedule( guided ) + for( unsigned idx_sol = 0; idx_sol < _setSize; idx_sol++ ) { + _optimizer->ReconstructMoments( _uSol[idx_sol], _alpha[idx_sol], _momentBasis ); + } } } diff --git a/src/datagenerator/datageneratorclassification.cpp b/src/datagenerator/datageneratorclassification.cpp index 666a599c..96dfbd5f 100644 --- a/src/datagenerator/datageneratorclassification.cpp +++ b/src/datagenerator/datageneratorclassification.cpp @@ -24,7 +24,6 @@ DataGeneratorClassification::DataGeneratorClassification( Config* settings ) : D _quadrature->ScalePointsAndWeights( _maxVelocity ); _optimizer->ScaleQuadWeights( _maxVelocity ); _weights = _quadrature->GetWeights(); - // std::cout << "sum of weights: " << _quadrature->SumUpWeights() << "\n"; _quadPointsSphere = _quadrature->GetPointsSphere(); _quadPoints = _quadrature->GetPoints(); @@ -99,12 +98,6 @@ Vector DataGeneratorClassification::ComputeMaxwellian( double rho, double u, dou Vector maxwellianAlpha = Vector( _nTotalEntries, 0.0 ); _optimizer->Solve( maxwellianAlpha, maxwellianMoment, _momentBasis ); - // std::cout << "Maxwellian Moment:\n"; - // std::cout << maxwellianMoment << std::endl; - // std::cout << moment0 << std::endl; - // std::cout << "Maxwellian Alpha:\n"; - // std::cout << maxwellianAlpha << std::endl; - maxwellianAlpha[1] *= 2; // For debugging, reconstruct the moments from Maxwellian alpha double entropyReconstruction = 0.0; @@ -113,8 +106,6 @@ Vector DataGeneratorClassification::ComputeMaxwellian( double rho, double u, dou entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( maxwellianAlpha, _momentBasis[idx_quad] ) ); maxwellianMoment += _momentBasis[idx_quad] * ( _weights[idx_quad] * entropyReconstruction ); } - // std::cout << "Reconstructed Maxwellian Moment:\n"; - // std::cout << maxwellianMoment << std::endl; return maxwellian; } @@ -124,6 +115,5 @@ double DataGeneratorClassification::ComputeKLDivergence( Vector& f1, Vector& f2 for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { sum += f1[idx_quad] * log( f1[idx_quad] / f2[idx_quad] ); } - // std::cout << "KL-Divergence:" << sum << std::endl; return sum; } diff --git a/src/datagenerator/datageneratorclassification2D.cpp b/src/datagenerator/datageneratorclassification2D.cpp index fc07f8f1..f463e711 100644 --- a/src/datagenerator/datageneratorclassification2D.cpp +++ b/src/datagenerator/datageneratorclassification2D.cpp @@ -250,21 +250,6 @@ void DataGeneratorClassification2D::PrintTrainingData() { for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { - // Test the moments of the kinetic density - // Vector u = Vector( _nTotalEntries, 0.0 ); - // double tmp = 0.0; - // double recT = 0.0; - // for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { - // u += _momentBasis[idx_quad] * _weights[idx_quad] * _kineticDensity[idx_set][idx_quad]; - // tmp += _weights[idx_quad] * _kineticDensity[idx_set][idx_quad]; - // recT += 0.5 * ( _momentBasis[idx_quad][1] * _momentBasis[idx_quad][1] + _momentBasis[idx_quad][2] * _momentBasis[idx_quad][2] ) * - // _weights[idx_quad] * _kineticDensity[idx_set][idx_quad]; - //} - // std::cout << "alpha_" << idx_set << " : " << _alpha[idx_set] << "\n"; - // std::cout << "u_" << idx_set << " : " << u << "\n"; - // std::cout << "rho" << idx_set << " : " << tmp << "\n"; - // std::cout << "T" << idx_set << " : " << recT << "\n"; - std::stringstream streamDensity; for( unsigned idx_quad = 0; idx_quad < _nq - 1; idx_quad++ ) { streamDensity << std::fixed << std::setprecision( 12 ) << _kineticDensity[idx_set][idx_quad] << ","; diff --git a/src/datagenerator/datageneratorregression.cpp b/src/datagenerator/datageneratorregression.cpp index a72f8728..d3e35e08 100644 --- a/src/datagenerator/datageneratorregression.cpp +++ b/src/datagenerator/datageneratorregression.cpp @@ -44,6 +44,7 @@ void DataGeneratorRegression::ComputeTrainingData() { VectorVector rot_uSol = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) ); VectorVector rot_alpha_comp = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) ); VectorVector rot_alpha = VectorVector( _setSize, Vector( _nTotalEntries, 0.0 ) ); + Vector alpha_norm_dummy( _setSize, 0 ); for( unsigned i = 0; i < _setSize; i++ ) { Vector u1{ _uSol[i][1], _uSol[i][2] }; @@ -70,7 +71,7 @@ void DataGeneratorRegression::ComputeTrainingData() { rot_alpha[i][3] = (float)( alpha2( 0, 1 ) ); rot_alpha[i][4] = (float)( alpha2( 1, 1 ) ); } - _optimizer->SolveMultiCell( rot_alpha_comp, rot_uSol, _momentBasis ); + _optimizer->SolveMultiCell( rot_alpha_comp, rot_uSol, _momentBasis, alpha_norm_dummy ); TextProcessingToolbox::PrintVectorVector( rot_alpha ); TextProcessingToolbox::PrintVectorVector( rot_alpha_comp ); @@ -84,74 +85,11 @@ void DataGeneratorRegression::ComputeTrainingData() { // ---- Check realizability --- // CheckRealizability(); - /* - std::vector timings( 100, 0.0 ); - for( unsigned j = 0; j < _setSize; j++ ) { - _uSol[j][0] = 1.0; - _uSol[j][1] = 0.99; - _uSol[j][2] = 0.99 * 0.99; - } - - for( unsigned i = 0; i < 100; i++ ) { - // Record start time - auto start = std::chrono::high_resolution_clock::now(); // --- compute alphas --- - _optimizer->SolveMultiCell( _alpha, _uSol, _momentBasis ); - - // Record end time - auto finish = std::chrono::high_resolution_clock::now(); - std::chrono::duration elapsed = finish - start; - timings[i] = elapsed.count(); - log->info( "| Elapsed time for solving the entropy minimization problem: " + std::to_string( elapsed.count() ) + " seconds. Iteration" + - std::to_string( i ) + "." ); - // reset training solution - for( unsigned j = 0; j < _setSize; j++ ) { - for( unsigned k = 0; k < _nTotalEntries; k++ ) { - _alpha[j][k] = 0.0; - } - } - } + Vector alpha_norm_dummy( _setSize, 0 ); - // Compute statistics - double mean = 0; - for( unsigned i = 0; i < 100; i++ ) { - mean += timings[i]; - } - mean /= 100; - double stdDev = 0; - for( unsigned i = 0; i < 100; i++ ) { - stdDev += ( timings[i] - mean ) * ( timings[i] - mean ); - } - stdDev /= 100; - stdDev = sqrt( stdDev ); - log->info( "| Mean timing: " + std::to_string( mean ) + " seconds. Standard deviation: " + std::to_string( stdDev ) + " seconds." ); - */ - - // some checks - // Vector u = { 1.0, 0.0, 0.5, 0.0 }; - // Vector alpha = { 1.0, 0.0, 0.5, 0.01 }; - //_optimizer->Solve( alpha, u, _momentBasis, 0 ); - // std::cout << "u=" << u << "\nalpha=" << alpha << "\n"; - // u = { 1.0, 0.0, 0.5, 0.01 }; - // alpha = { 1.0, 0.0, 0.5, 0.01 }; - //_optimizer->Solve( alpha, u, _momentBasis, 0 ); - // std::cout << "u=" << u << "\nalpha=" << alpha << "\n"; - // u = { 1.0, 0.0, 0.5, 0.05 }; - // alpha = { 1.0, 0.0, 0.5, 0.01 }; - //_optimizer->Solve( alpha, u, _momentBasis, 0 ); - // std::cout << "u=" << u << "\nalpha=" << alpha << "\n"; - // u = { 1.0, 0.0, 0.5, 0.1 }; - // alpha = { 1.0, 0.0, 0.5, 0.01 }; - //_optimizer->Solve( alpha, u, _momentBasis, 0 ); - // std::cout << "u=" << u << "\nalpha=" << alpha << "\n"; - // u = { 1.0, 0.0, 0.5, 0.2 }; - // alpha = { 1.0, 0.0, 0.5, 0.01 }; - //_optimizer->Solve( alpha, u, _momentBasis, 0 ); - // std::cout << "u=" << u << "\nalpha=" << alpha << "\n"; - - // --- compute alphas --- - _optimizer->SolveMultiCell( _alpha, _uSol, _momentBasis ); + _optimizer->SolveMultiCell( _alpha, _uSol, _momentBasis, alpha_norm_dummy ); // --- Postprocessing if( _settings->GetRealizabilityReconstruction() ) { log->info( "| Compute realizable problems." ); @@ -195,27 +133,48 @@ void DataGeneratorRegression::ComputeEntropyH_dual() { } void DataGeneratorRegression::ComputeEntropyH_primal() { + if( _reducedSampling ) { + VectorVector momentsRed = VectorVector( _nq, Vector( _nTotalEntries - 1, 0.0 ) ); + + for( unsigned idx_nq = 0; idx_nq < _nq; idx_nq++ ) { // copy (reduced) moments + for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { + momentsRed[idx_nq][idx_sys - 1] = _momentBasis[idx_nq][idx_sys]; + } + } #pragma omp parallel for schedule( guided ) - for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { - _hEntropy[idx_set] = -1 * _optimizer->ComputeObjFunc( _alpha[idx_set], _uSol[idx_set], _momentBasis ); + for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { + Vector uSolRed( _nTotalEntries - 1, 0.0 ); + Vector alphaRed( _nTotalEntries - 1, 0.0 ); + for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) { + alphaRed[idx_sys - 1] = _alpha[idx_set][idx_sys]; + uSolRed[idx_sys - 1] = _uSol[idx_set][idx_sys]; + } + _hEntropy[idx_set] = -1 * _optimizer->ComputeObjFunc( alphaRed, uSolRed, momentsRed ); + } } + else { +#pragma omp parallel for schedule( guided ) + for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { + _hEntropy[idx_set] = -1 * _optimizer->ComputeObjFunc( _alpha[idx_set], _uSol[idx_set], _momentBasis ); + } - for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { - if( isnan( _hEntropy[idx_set] ) ) { - std::string msgU = "("; - std::string msgAlpha = "("; - for( unsigned idx_basis = 0; idx_basis < _nTotalEntries - 1; idx_basis++ ) { - msgU += std::to_string( _uSol[idx_set][idx_basis] ) + "|"; - msgAlpha += std::to_string( _alpha[idx_set][idx_basis] ) + "|"; + for( unsigned idx_set = 0; idx_set < _setSize; idx_set++ ) { + if( isnan( _hEntropy[idx_set] ) ) { + std::string msgU = "("; + std::string msgAlpha = "("; + for( unsigned idx_basis = 0; idx_basis < _nTotalEntries - 1; idx_basis++ ) { + msgU += std::to_string( _uSol[idx_set][idx_basis] ) + "|"; + msgAlpha += std::to_string( _alpha[idx_set][idx_basis] ) + "|"; + } + msgU += std::to_string( _uSol[idx_set][_nTotalEntries - 1] ) + ")"; + msgAlpha += std::to_string( _uSol[idx_set][_nTotalEntries - 1] ) + ")"; + + ErrorMessages::Error( "Value for h is NaN. This can happen near the boundary of the realizable set.\nu= " + msgU + + "\nalpha= " + msgAlpha + + "\nPlease adjust the Options " + "REALIZABLE_SET_EPSILON_U0 and REALIZABLE_SET_EPSILON_U1.", + CURRENT_FUNCTION ); } - msgU += std::to_string( _uSol[idx_set][_nTotalEntries - 1] ) + ")"; - msgAlpha += std::to_string( _uSol[idx_set][_nTotalEntries - 1] ) + ")"; - - ErrorMessages::Error( "Value for h is NaN. This can happen near the boundary of the realizable set.\nu= " + msgU + - "\nalpha= " + msgAlpha + - "\nPlease adjust the Options " - "REALIZABLE_SET_EPSILON_U0 and REALIZABLE_SET_EPSILON_U1.", - CURRENT_FUNCTION ); } } } diff --git a/src/datagenerator/datageneratorregression3D.cpp b/src/datagenerator/datageneratorregression3D.cpp index ad69fe47..c3e936d4 100644 --- a/src/datagenerator/datageneratorregression3D.cpp +++ b/src/datagenerator/datageneratorregression3D.cpp @@ -137,9 +137,7 @@ void DataGeneratorRegression3D::CheckRealizability() { CURRENT_FUNCTION ); } if( normU1 / _uSol[idx_set][0] <= 0 /*+ 0.5 * epsilon*/ ) { - // std::cout << "_uSol" << _uSol[idx_set][1] << " | " << _uSol[idx_set][2] << " | " << _uSol[idx_set][3] << " \n"; - // std::cout << "normU1 / _uSol[" << idx_set << "][0]: " << normU1 / _uSol[idx_set][0] << "\n"; - // std::cout << "normU1: " << normU1 << " | _uSol[idx_set][0] " << _uSol[idx_set][0] << "\n"; + ErrorMessages::Error( "Moment to close to boundary of realizable set [code 2].\nBoundary ratio: " + std::to_string( normU1 / _uSol[idx_set][0] ), CURRENT_FUNCTION ); diff --git a/src/optimizers/neuralnetworkoptimizer.cpp b/src/optimizers/neuralnetworkoptimizer.cpp index b89d4ec6..6e121103 100644 --- a/src/optimizers/neuralnetworkoptimizer.cpp +++ b/src/optimizers/neuralnetworkoptimizer.cpp @@ -130,7 +130,7 @@ Vector NeuralNetworkOptimizer::RotateM1( Vector& vec, Matrix& R ) { return R * v Matrix NeuralNetworkOptimizer::RotateM2( Matrix& vec, Matrix& R, Matrix& Rt ) { return R * vec * Rt; } -void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& /*moments*/ ) { +void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& u, const VectorVector& /*moments*/, Vector& alpha_norms ) { unsigned servingSize = _settings->GetNCells(); Matrix rot180{ { -1.0, 0.0 }, { 0.0, -1.0 } }; @@ -156,7 +156,8 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } else { // Rotate everything to x-Axis of first moment tensor - //#pragma omp parallel for + //std::cout << "rotation_active\n"; +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _settings->GetNCells(); idx_cell++ ) { Vector u1{ u[idx_cell][1], u[idx_cell][2] }; Matrix u2{ { u[idx_cell][3], u[idx_cell][4] }, { u[idx_cell][4], u[idx_cell][5] } }; @@ -253,6 +254,7 @@ void NeuralNetworkOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& } alphaRed[idx_sys] = ( alphaRed[idx_sys] + alphaRedMirror[idx_sys] ) / 2; // average (and store in alphaRed) } + alpha_norms[idx_cell] = norm( alphaRed ) * norm( alphaRed ); // Rotate Back Vector alpha1{ alphaRed[0], alphaRed[1] }; diff --git a/src/optimizers/newtonoptimizer.cpp b/src/optimizers/newtonoptimizer.cpp index d9adaf12..54d57fa6 100644 --- a/src/optimizers/newtonoptimizer.cpp +++ b/src/optimizers/newtonoptimizer.cpp @@ -67,9 +67,10 @@ void NewtonOptimizer::ComputeHessian( Vector& alpha, const VectorVector& moments } } -void NewtonOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& sol, const VectorVector& moments ) { +void NewtonOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& sol, const VectorVector& moments, Vector& alpha_norms ) { unsigned nCells = alpha.size(); + unsigned nSys = alpha[0].size(); // if we have quadratic entropy, then alpha = u; if( _settings->GetEntropyName() == QUADRATIC && _settings->GetNewtonFastMode() ) { @@ -82,6 +83,8 @@ void NewtonOptimizer::SolveMultiCell( VectorVector& alpha, VectorVector& sol, co #pragma omp parallel for schedule( guided ) for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) { Solve( alpha[idx_cell], sol[idx_cell], moments, idx_cell ); + for ( unsigned idx_sys = 1; idx_sys < nSys; idx_sys++ ) + alpha_norms[idx_cell] *= alpha[idx_cell][idx_sys]*alpha[idx_cell][idx_sys]; } } @@ -158,7 +161,7 @@ void NewtonOptimizer::Solve( Vector& alpha, Vector& sol, const VectorVector& mom alphaNew = alpha - stepSize * _delta * H * grad; ComputeGradient( alphaNew, sol, moments, dalphaNew ); - // Check if FONC is fullfilled + // Check if FONC is locally fullfilled if( norm( dalphaNew ) < _epsilon ) { alpha = alphaNew; return; diff --git a/src/optimizers/partregularizednewtonoptimizer.cpp b/src/optimizers/partregularizednewtonoptimizer.cpp index 38e3d446..f170d661 100644 --- a/src/optimizers/partregularizednewtonoptimizer.cpp +++ b/src/optimizers/partregularizednewtonoptimizer.cpp @@ -8,13 +8,10 @@ #include "optimizers/partregularizednewtonoptimizer.hpp" #include "common/config.hpp" #include "entropies/entropybase.hpp" -#include "optimizers/newtonoptimizer.hpp" #include "quadratures/quadraturebase.hpp" #include "toolboxes/errormessages.hpp" #include "toolboxes/textprocessingtoolbox.hpp" -#include - PartRegularizedNewtonOptimizer::PartRegularizedNewtonOptimizer( Config* settings ) : NewtonOptimizer( settings ) { _gamma = settings->GetRegularizerGamma(); // Regularization parameter (to be included in settings) } diff --git a/src/optimizers/reducednewtonoptimizer.cpp b/src/optimizers/reducednewtonoptimizer.cpp new file mode 100644 index 00000000..577c4ab4 --- /dev/null +++ b/src/optimizers/reducednewtonoptimizer.cpp @@ -0,0 +1,111 @@ +/*! + * @file reducednewtonoptimizer.cpp + * @brief class for solving the minimal entropy optimization problem using a newton optimizer with line search. + * @author S. Schotthöfer + */ + +#include "optimizers/reducednewtonoptimizer.hpp" +#include "common/config.hpp" +#include "entropies/entropybase.hpp" +#include "quadratures/quadraturebase.hpp" +#include "toolboxes/errormessages.hpp" +#include "toolboxes/textprocessingtoolbox.hpp" + +ReducedNewtonOptimizer::ReducedNewtonOptimizer( Config* settings ) : NewtonOptimizer( settings ) {} + +ReducedNewtonOptimizer::~ReducedNewtonOptimizer() { delete _quadrature; } + +double ReducedNewtonOptimizer::ComputeObjFunc( Vector& alpha, Vector& sol, const VectorVector& moments ) { + double result = 0.0; + + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + result += _entropy->EntropyDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad]; + } + // log + result = log( result ); + // linear term + result -= dot( alpha, sol ); + return result + 1.0; +} + +void ReducedNewtonOptimizer::ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) { + + // Reset Vector + for( unsigned idx_sys = 0; idx_sys < grad.size(); idx_sys++ ) { + grad[idx_sys] = 0.0; + } + + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + grad += moments[idx_quad] * ( _entropy->EntropyPrimeDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] ); + } + + double result = 0.0; + + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + result += _entropy->EntropyDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad]; + } + + grad /= result; // normalization + + // linear part + grad -= sol; +} + +void ReducedNewtonOptimizer::ComputeHessian( Vector& alpha, const VectorVector& moments, Matrix& hessian ) { + // Reset Matrix + unsigned nSize = alpha.size(); + + for( unsigned idx_Row = 0; idx_Row < nSize; idx_Row++ ) { + for( unsigned idx_Col = 0; idx_Col < nSize; idx_Col++ ) { + hessian( idx_Row, idx_Col ) = 0.0; + // if( idx_Col == idx_Row ) hessian( idx_Row, idx_Col ) = 1.0; + } + } + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + hessian += + outer( moments[idx_quad], moments[idx_quad] ) * ( _entropy->EntropyHessianDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] ); + } + + double result = 0.0; + + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + result += _entropy->EntropyDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad]; + } + + hessian /= result; // normalization + + // second part + Vector grad( nSize, 0.0 ); + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + grad += moments[idx_quad] * ( _entropy->EntropyPrimeDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] ); + } + + hessian += -1 * outer( grad, grad ) / ( result * result ); +} + +void ReducedNewtonOptimizer::ReconstructMoments( Vector& sol, const Vector& alpha, const VectorVector& moments ) { + double entropyReconstruction = 0.0; + for( unsigned idx_sys = 0; idx_sys < sol.size(); idx_sys++ ) { + sol[idx_sys] = 0.0; + } + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + // Make entropyReconstruction a member vector, s.t. it does not have to be re-evaluated in ConstructFlux + entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( alpha, moments[idx_quad] ) ); + sol += moments[idx_quad] * ( _weights[idx_quad] * entropyReconstruction ); + } + + double result = 0.0; + + // Integrate + for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { + result += _entropy->EntropyDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad]; + } + + sol /= result; // normalization +} diff --git a/src/optimizers/reducedpartregularizednewtonoptimizer.cpp b/src/optimizers/reducedpartregularizednewtonoptimizer.cpp new file mode 100644 index 00000000..04c3a0fb --- /dev/null +++ b/src/optimizers/reducedpartregularizednewtonoptimizer.cpp @@ -0,0 +1,48 @@ +/*! + * @file reducedpartregularizednewtonoptimizer.cpp + * @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 + */ + +#include "optimizers/reducedpartregularizednewtonoptimizer.hpp" +#include "common/config.hpp" +#include "entropies/entropybase.hpp" +#include "quadratures/quadraturebase.hpp" +#include "toolboxes/errormessages.hpp" +#include "toolboxes/textprocessingtoolbox.hpp" + +ReducedPartRegularizedNewtonOptimizer::ReducedPartRegularizedNewtonOptimizer( Config* settings ) : ReducedNewtonOptimizer( settings ) { + _gamma = settings->GetRegularizerGamma(); // Regularization parameter (to be included in settings) +} + +ReducedPartRegularizedNewtonOptimizer::~ReducedPartRegularizedNewtonOptimizer() {} + +double ReducedPartRegularizedNewtonOptimizer::ComputeObjFunc( Vector& alpha, Vector& sol, const VectorVector& moments ) { + double result = ReducedNewtonOptimizer::ComputeObjFunc( alpha, sol, moments ); // Calls non regularized objective function + for( unsigned idx_sys = 0; idx_sys < alpha.size(); idx_sys++ ) { + result += 0.5 * _gamma * alpha[idx_sys] * alpha[idx_sys]; // Add regularizer norm(_alpha^r)^2 + } + return result; +} + +void ReducedPartRegularizedNewtonOptimizer::ComputeGradient( Vector& alpha, Vector& sol, const VectorVector& moments, Vector& grad ) { + ReducedNewtonOptimizer::ComputeGradient( alpha, sol, moments, grad ); // compute unregularized gradients + for( unsigned idx_sys = 0; idx_sys < alpha.size(); idx_sys++ ) { + grad[idx_sys] += _gamma * alpha[idx_sys]; // Add regularizer _alpha^r + } +} + +void ReducedPartRegularizedNewtonOptimizer::ComputeHessian( Vector& alpha, const VectorVector& moments, Matrix& hessian ) { + ReducedNewtonOptimizer::ComputeHessian( alpha, moments, hessian ); // compute unregularized hessian) + for( unsigned idx_sys = 0; idx_sys < alpha.size(); idx_sys++ ) { + hessian( idx_sys, idx_sys ) += _gamma; // Add block identity matrix with regularizer + } +} + +void ReducedPartRegularizedNewtonOptimizer::ReconstructMoments( Vector& sol, const Vector& alpha, const VectorVector& moments ) { + ReducedNewtonOptimizer::ReconstructMoments( sol, alpha, moments ); + for( unsigned idx_sys = 0; idx_sys < alpha.size(); idx_sys++ ) { + sol[idx_sys] += _gamma * alpha[idx_sys]; // Add regularizer _alpha^r + } +} diff --git a/src/problems/starmapvalidation.cpp b/src/problems/starmapvalidation.cpp index 7ee34dda..bc50a749 100644 --- a/src/problems/starmapvalidation.cpp +++ b/src/problems/starmapvalidation.cpp @@ -88,7 +88,6 @@ StarMapValidation_Moment::StarMapValidation_Moment( Config* settings, Mesh* mesh StarMapValidation_Moment::~StarMapValidation_Moment() {} VectorVector StarMapValidation_Moment::SetupIC() { - SphericalBase* tempBase = SphericalBase::Create( _settings ); unsigned ntotalEquations = tempBase->GetBasisSize(); delete tempBase; @@ -103,6 +102,7 @@ VectorVector StarMapValidation_Moment::SetupIC() { x = cellMidpoints[idx_cell][0]; y = cellMidpoints[idx_cell][1]; f = NormPDF( x, pos_beam[0], stddev ) * NormPDF( y, pos_beam[1], stddev ); + initialSolution[idx_cell][0] = f * StarMAPmoments[0]; for( unsigned idx_sys = 1; idx_sys < ntotalEquations; idx_sys++ ) { initialSolution[idx_cell][idx_sys] = f * StarMAPmoments[idx_sys]; // must be VectorVector diff --git a/src/solvers/csdmnsolver.cpp b/src/solvers/csdmnsolver.cpp index 563ed465..01025aba 100644 --- a/src/solvers/csdmnsolver.cpp +++ b/src/solvers/csdmnsolver.cpp @@ -71,11 +71,12 @@ CSDMNSolver::~CSDMNSolver() {} void CSDMNSolver::IterPreprocessing( unsigned idx_iter ) { // ------- Entropy closure Step ---------------- + Vector alpha_norm_dummy( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE - _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis ); + _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis, alpha_norm_dummy ); // ------- Solution reconstruction step ---- - //#pragma omp parallel for +#pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { // compute the kinetic density at all grid cells diff --git a/src/solvers/mnsolver.cpp b/src/solvers/mnsolver.cpp index 9c502185..e8882c1e 100644 --- a/src/solvers/mnsolver.cpp +++ b/src/solvers/mnsolver.cpp @@ -110,7 +110,9 @@ void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) { void MNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) { // ------- Entropy closure Step ---------------- - _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis ); // parallel + Vector alpha_norm_dummy( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE + + _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis, alpha_norm_dummy ); // parallel // ------- Solution reconstruction step ---- #pragma omp parallel for diff --git a/src/solvers/mnsolver_normalized.cpp b/src/solvers/mnsolver_normalized.cpp index 2743c415..47559aed 100644 --- a/src/solvers/mnsolver_normalized.cpp +++ b/src/solvers/mnsolver_normalized.cpp @@ -30,15 +30,20 @@ void MNSolverNormalized::IterPreprocessing( unsigned /*idx_pseudotime*/ ) { } // TextProcessingToolbox::PrintVectorVector( _sol ); - _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis ); + Vector alpha_norm_per_cell( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE + + _optimizer->SolveMultiCell( _alpha, _sol, _momentBasis, alpha_norm_per_cell ); // ------- Solution reconstruction step ---- #pragma omp parallel for for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { + + alpha_norm_per_cell[idx_cell] *= _momentBasis[0][0] * 0.5 * _settings->GetRegularizerGamma(); // is constant + // std::cout << alpha_norm << "|" << _momentBasis[0][0] << "\n"; for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) { // compute the kinetic density at all grid cells _kineticDensity[idx_cell][idx_quad] = - _u0[idx_cell] * _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _momentBasis[idx_quad] ) ); + _u0[idx_cell] * _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _momentBasis[idx_quad] ) - alpha_norm_per_cell[idx_cell] ); } if( _settings->GetRealizabilityReconstruction() ) ComputeRealizableSolution( idx_cell ); _sol[idx_cell] *= _u0[idx_cell]; diff --git a/src/solvers/pnsolver.cpp b/src/solvers/pnsolver.cpp index 73805fdb..55b7535a 100644 --- a/src/solvers/pnsolver.cpp +++ b/src/solvers/pnsolver.cpp @@ -226,7 +226,7 @@ void PNSolver::FVMUpdate( unsigned idx_energy ) { - _dE * _sol[idx_cell][idx_sys] * ( _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] /* scattering influence */ + _sigmaT[idx_energy][idx_cell] ); - /* total xs influence */ // Vorzeichenfehler! + /* total xs influence */ } // Source Term _solNew[idx_cell] += _dE * _Q[0][idx_cell]; diff --git a/tools/CI/docker_run_interactive_specify_cores.sh b/tools/CI/docker_run_interactive_specify_cores.sh index 94c11cb5..bea10910 100644 --- a/tools/CI/docker_run_interactive_specify_cores.sh +++ b/tools/CI/docker_run_interactive_specify_cores.sh @@ -1 +1 @@ -docker run -i -t --rm --cpuset-cpus=0-8 -v $(pwd)/../..:/mnt kitrt/test_ml:latest /bin/bash +docker run -i -t --rm --cpuset-cpus=0-12 -v $(pwd)/../..:/mnt kitrt/test_ml:latest /bin/bash