CSMMLab / KiT-RT / #95

Committed 28 May 2025 02:06AM UTC coverage: 46.655% (-11.1%) from 57.728%

Build # #95

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

Release 2025 (#44)

* New neural models (#30)

* extend kitrt script

* added new neural models

* extend to m3 models

* added M3_2D models

* added M2 and M4 models

* configure ml optimizer for high order models

* added configuration for high order models

* add option for rotation postprcessing in neural networks. remove unneccessary python scripts

* started rotational symmetric postprocessing

* added rotational invariance postprocessing for m2 and m1 models

* fix post merge bugs

* created hohlraum mesh

* add hohlraum tes case

* add hohlraum test case

* add hohlraum cfg filees

* fixed hohlraum testcase

* add hohlraum cfg files

* changed hohlraum cfg

* changed hohlraum testcase to isotropic inflow source boundary condition

* added ghost cell bonudary for hohlraum testcase

* update readme and linesource mesh creator

* added proper scaling for linesource reference solution

* regularized newton debugging

* Data generator with reduced objective functional (#33)

* added reduced optimizer for sampling

* remove old debugging comments

* mesh acceleration (#38)

* added new ansatz (#36)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* Mesh acc (#41)

* mesh acceleration

* added floor value for starmap moment methods

* enable accelleration

* delete minimum value starmap

* Update README.md

* Spherical harmonics nn (#40)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* added sph... (continued)

Run Details

767 of 3019 new or added lines in 51 files covered. (25.41%)

131 existing lines in 8 files now uncovered.

4422 of 9478 relevant lines covered (46.66%)

57163.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

85.29

/src/solvers/mnsolver.cpp

#include "solvers/mnsolver.hpp"
#include "common/config.hpp"
#include "common/io.hpp"
#include "common/mesh.hpp"
#include "entropies/entropybase.hpp"
#include "fluxes/numericalflux.hpp"
#include "optimizers/optimizerbase.hpp"
#include "problems/problembase.hpp"
#include "quadratures/quadraturebase.hpp"
#include "toolboxes/errormessages.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"
#include "velocitybasis/sphericalbase.hpp"

// externals
#include "spdlog/spdlog.h"

MNSolver::MNSolver( Config* settings ) : SolverBase( settings ) {

    _polyDegreeBasis = settings->GetMaxMomentDegree();
    _basis           = SphericalBase::Create( _settings );
    _nSystem         = _basis->GetBasisSize();

    // build quadrature object and store quadrature points and weights
    _quadPoints = _quadrature->GetPoints();
    _weights    = _quadrature->GetWeights();
    //_nq               = _quadrature->GetNq();
    _quadPointsSphere = _quadrature->GetPointsSphere();
    //_settings->SetNQuadPoints( _nq );

    // Initialize Entropy
    _entropy = EntropyBase::Create( _settings );

    // Initialize Optimizer
    _optimizer = OptimizerBase::Create( _settings );

    // Initialize lagrange Multiplier
    _alpha = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );

    // Initialize kinetic density at grid cells
    _kineticDensity = VectorVector( _nCells, Vector( _nq, 0.0 ) );

    // Limiter variables
    _solDx   = VectorVector( _nCells, Vector( _nq, 0.0 ) );
    _solDy   = VectorVector( _nCells, Vector( _nq, 0.0 ) );
    _limiter = VectorVector( _nCells, Vector( _nq, 0.0 ) );

    // Initialize and Pre-Compute Moments at quadrature points
    _momentBasis = VectorVector( _nq, Vector( _nSystem, 0.0 ) );
    ComputeMoments();

    // Initialize Scatter Matrix --
    _scatterMatDiag = Vector( _nSystem, 0.0 );
    ComputeScatterMatrix();
}

MNSolver::~MNSolver() {
    delete _entropy;
    delete _optimizer;
    delete _basis;
}

void MNSolver::ComputeScatterMatrix() {

    // --- Isotropic ---
    if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
        _scatterMatDiag[0] = 1.0;
        for( unsigned idx_diag = 1; idx_diag < _nSystem; idx_diag++ ) {
            _scatterMatDiag[idx_diag] = 0.0;    // SPHERICAL_HARMONICS are orthogonal basis
        }
    }
    else {    // SPHERICAL_MONOMIALS
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
            _scatterMatDiag[idx_sys] = 0.0;
            for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
                _scatterMatDiag[idx_sys] += _momentBasis[idx_quad][idx_sys] * _weights[idx_quad];
            }
            if( _settings->GetDim() == 1 ) {
                _scatterMatDiag[idx_sys] /= 2;
            }
            else if( _settings->GetDim() == 2 ) {
                _scatterMatDiag[idx_sys] /= M_PI;
            }
            else {    // 3D
                _scatterMatDiag[idx_sys] /= 4.0 * M_PI;
            }
        }
    }
}

void MNSolver::ComputeMoments() {
    double my, phi;

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
        my                     = _quadPointsSphere[idx_quad][0];
        phi                    = _quadPointsSphere[idx_quad][1];
        _momentBasis[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
    }
}

void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
    // double entropyReconstruction = 0.0;

    for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {    // reset solution
        _sol[idx_cell][idx_sys] = 0.0;
    }
    _optimizer->ReconstructMoments( _sol[idx_cell], _alpha[idx_cell], _momentBasis );
}

void MNSolver::IterPreprocessing( unsigned /*idx_pseudotime*/ ) {

    // ------- Entropy closure Step ----------------
    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
    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
            _kineticDensity[idx_cell][idx_quad] = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _momentBasis[idx_quad] ) );
        }
        if( _settings->GetRealizabilityReconstruction() ) ComputeRealizableSolution( idx_cell );
    }
    // ------ Compute slope limiters and cell gradients ---
    if( _reconsOrder > 1 ) {
        _mesh->ComputeSlopes( _nq, _solDx, _solDy, _kineticDensity );               // parallel
        _mesh->ComputeLimiter( _nq, _solDx, _solDy, _kineticDensity, _limiter );    // parallel
    }
}

void MNSolver::ComputeScalarFlux() {
    double firstMomentScaleFactor = 4 * M_PI;
    if( _settings->GetProblemName() == PROBLEM_Aircavity1D || _settings->GetProblemName() == PROBLEM_Linesource1D ||
        _settings->GetProblemName() == PROBLEM_Checkerboard1D ) {
        firstMomentScaleFactor = 2.0;
    }
    if( _settings->GetDim() == 2 ) {
        firstMomentScaleFactor = M_PI;
    }
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        _scalarFluxNew[idx_cell] = _sol[idx_cell][0] * firstMomentScaleFactor;
    }
}

void MNSolver::FluxUpdate() {
    // Loop over all spatial cells
    if( _settings->GetProblemName() == PROBLEM_Aircavity1D || _settings->GetProblemName() == PROBLEM_Linesource1D ||
        _settings->GetProblemName() == PROBLEM_Checkerboard1D ) {
        FluxUpdatePseudo1D();
    }
    else {
        FluxUpdatePseudo2D();
    }
}

void MNSolver::FluxUpdatePseudo1D() {
// Loop over the grid cells
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        // Dirichlet Boundaries stayd
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        double solL, solR, kineticFlux;
        Vector flux( _nSystem, 0.0 );
        for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
            kineticFlux = 0.0;    // reset temorary flux

            for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); idx_nbr++ ) {
                if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
                    // Boundary cells are first order and mirror ghost cells
                    solL = _kineticDensity[idx_cell][idx_quad];
                    solR = solL;
                }
                else {
                    // interior cell
                    unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr];    // global idx of neighbor cell
                    if( _reconsOrder == 1 ) {
                        solL = _kineticDensity[idx_cell][idx_quad];
                        solR = _kineticDensity[_neighbors[idx_cell][idx_nbr]][idx_quad];
                    }
                    else if( _reconsOrder == 2 ) {
                        solL = _kineticDensity[idx_cell][idx_quad] +
                               _limiter[idx_cell][idx_quad] *
                                   ( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
                                     _solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
                        solR = _kineticDensity[nbr_glob][idx_quad] +
                               _limiter[nbr_glob][idx_quad] *
                                   ( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
                                     _solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
                    }
                    else {
                        ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION );
                    }
                }
                // Kinetic flux
                kineticFlux += _g->Flux1D( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
            }
            // Solution flux
            flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
        }
        _solNew[idx_cell] = flux;
    }
}

void MNSolver::FluxUpdatePseudo2D() {
// Loop over the grid cells
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        // Dirichlet Boundaries stayd
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        double solL, solR, kineticFlux;
        Vector flux( _nSystem, 0.0 );
        for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
            kineticFlux = 0.0;    // reset temorary flux
            for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); idx_nbr++ ) {
                if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
                    // Boundary cells are first order and mirror ghost cells
                    solL = _kineticDensity[idx_cell][idx_quad];
                    solR = solL;
                }
                else {
                    // interior cell
                    unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr];    // global idx of neighbor cell
                    if( _reconsOrder == 1 ) {
                        solL = _kineticDensity[idx_cell][idx_quad];
                        solR = _kineticDensity[_neighbors[idx_cell][idx_nbr]][idx_quad];
                    }
                    else if( _reconsOrder == 2 ) {
                        solL = _kineticDensity[idx_cell][idx_quad] +
                               _limiter[idx_cell][idx_quad] *
                                   ( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
                                     _solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
                        solR = _kineticDensity[nbr_glob][idx_quad] +
                               _limiter[nbr_glob][idx_quad] *
                                   ( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
                                     _solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
                    }
                    else {
                        ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION );
                    }
                }
                // Kinetic flux
                kineticFlux += _g->Flux( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
            }
            // Solution flux
            flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
        }
        _solNew[idx_cell] = flux;
    }
}

void MNSolver::FVMUpdate( unsigned idx_iter ) {
    if( _Q.size() == 1u && _sigmaT.size() == 1u && _sigmaS.size() == 1u ) {    // Physics constant in time
// Loop over the grid cells
#pragma omp parallel for
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
            // Dirichlet Boundaries stay
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
            // Flux update
            for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
                _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys]                                                     /* old solution */
                                             - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys]                   /* cell averaged flux */
                                             + _dT * _sigmaS[0][idx_cell] * _sol[idx_cell][0] * _scatterMatDiag[idx_sys] /* scattering gain */
                                             - _dT * ( _sigmaT[0][idx_cell] ) * _sol[idx_cell][idx_sys]; /* scattering and absorbtion loss */
            }
            // Source Term
            _solNew[idx_cell] += _dT * _Q[0][idx_cell];
        }
    }
    else {
        // Loop over the grid cells
#pragma omp parallel for
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
            // Dirichlet Boundaries stay
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
            // Flux update
            for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
                _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys]                                   /* old solution */
                                             - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
                                             + _dT * _sigmaS[idx_iter][idx_cell] * _sol[idx_cell][0] * _scatterMatDiag[idx_sys] /* scattering gain */
                                             - _dT * ( _sigmaT[idx_iter][idx_cell] ) * _sol[idx_cell][idx_sys]; /* scattering and absorbtion loss */
            }
            // Source Term
            _solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell];
        }
    }
}

void MNSolver::PrepareVolumeOutput() {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
    _outputFieldNames.resize( nGroups );
    _outputFields.resize( nGroups );
    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
#pragma omp parallel for
    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
        // Prepare all Output Fields per group
        // Different procedure, depending on the Group...
        switch( _settings->GetVolumeOutput()[idx_group] ) {
            case MINIMAL:
                // Currently only one entry ==> rad flux
                _outputFields[idx_group].resize( 1 );
                _outputFieldNames[idx_group].resize( 1 );
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
                break;
            case MOMENTS:
                // As many entries as there are moments in the system
                _outputFields[idx_group].resize( _nSystem );
                _outputFieldNames[idx_group].resize( _nSystem );
                if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                    if( _settings->GetDim() == 3 ) {
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {    // 3D
                            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                                _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                                _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                    std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                            }
                        }
                    }
                    else if( _settings->GetDim() == 2 ) {
                        unsigned count = 0;
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                                if( ( idx_l + idx_k ) % 2 == 0 ) {
                                    _outputFields[idx_group][count].resize( _nCells );
                                    _outputFieldNames[idx_group][count] =
                                        std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                                    count++;
                                }
                            }
                        }
                    }
                    else if( _settings->GetDim() == 1 ) {
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                            _outputFields[idx_group][idx_l].resize( _nCells );
                            _outputFieldNames[idx_group][idx_l] = std::string( "u_" + std::to_string( idx_l ) + "^0" );
                        }
                    }
                }
                else {
                    for( unsigned idx_l = 0; idx_l <= _polyDegreeBasis; idx_l++ ) {
                        unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
                        for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
                            _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
                    }
                }
                break;
            case DUAL_MOMENTS:
                // As many entries as there are moments in the system
                _outputFields[idx_group].resize( _nSystem );
                _outputFieldNames[idx_group].resize( _nSystem );
                if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                    if( _settings->GetDim() == 3 ) {
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                                _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                                _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                    std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                            }
                        }
                    }
                    else if( _settings->GetDim() == 2 ) {
                        unsigned count = 0;
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                                if( ( idx_l + idx_k ) % 2 == 0 ) {
                                    _outputFields[idx_group][count].resize( _nCells );
                                    _outputFieldNames[idx_group][count] =
                                        std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                                    count++;
                                }
                            }
                        }
                    }
                    else if( _settings->GetDim() == 1 ) {
                        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                            _outputFields[idx_group][idx_l].resize( _nCells );
                            _outputFieldNames[idx_group][idx_l] = std::string( "alpha_" + std::to_string( idx_l ) + "^0" );
                        }
                    }
                }
                else {    // SPHERICAL_MONOMIALS
                    for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                        unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
                        for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
                            _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
                    }
                }
                break;
            case ANALYTIC:
                // one entry per cell
                _outputFields[idx_group].resize( 1 );
                _outputFieldNames[idx_group].resize( 1 );
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = std::string( "analytic radiation flux density" );
                break;
            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

void MNSolver::WriteVolumeOutput( unsigned idx_iter ) {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
    // Check if volume output fields are written to file this iteration
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
        ( idx_iter == _nEnergies - 1 ) /* need sol at last iteration */ ) {
#pragma omp parallel for
        for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
            switch( _settings->GetVolumeOutput()[idx_group] ) {
                case MINIMAL:
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        _outputFields[idx_group][0][idx_cell] = _scalarFluxNew[idx_cell];
                    }
                    break;
                case MOMENTS:
                    for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
                        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                            _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
                        }
                    }
                    break;
                case DUAL_MOMENTS:
                    for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
                        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
                            _outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];
                        }
                    }
                    break;
                case ANALYTIC:
                    // Compute total "mass" of the system ==> to check conservation properties
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        double time                           = idx_iter * _dT;
                        _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
                            _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_iter][idx_cell] );
                    }
                    break;
                default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
            }
        }
    }
}

1	#include "solvers/mnsolver.hpp"
2	#include "common/config.hpp"
3	#include "common/io.hpp"
4	#include "common/mesh.hpp"
5	#include "entropies/entropybase.hpp"
6	#include "fluxes/numericalflux.hpp"
7	#include "optimizers/optimizerbase.hpp"
8	#include "problems/problembase.hpp"
9	#include "quadratures/quadraturebase.hpp"
10	#include "toolboxes/errormessages.hpp"
11	#include "toolboxes/textprocessingtoolbox.hpp"
12	#include "velocitybasis/sphericalbase.hpp"
13
14	// externals
15	#include "spdlog/spdlog.h"
16
17	MNSolver::MNSolver( Config* settings ) : SolverBase( settings ) {	8✔
18
19	_polyDegreeBasis = settings->GetMaxMomentDegree();	8✔
20	_basis = SphericalBase::Create( _settings );	8✔
21	_nSystem = _basis->GetBasisSize();	8✔
22
23	// build quadrature object and store quadrature points and weights
24	_quadPoints = _quadrature->GetPoints();	8✔
25	_weights = _quadrature->GetWeights();	8✔
26	//_nq = _quadrature->GetNq();
27	_quadPointsSphere = _quadrature->GetPointsSphere();	8✔
28	//_settings->SetNQuadPoints( _nq );
29
30	// Initialize Entropy
31	_entropy = EntropyBase::Create( _settings );	8✔
32
33	// Initialize Optimizer
34	_optimizer = OptimizerBase::Create( _settings );	8✔
35
36	// Initialize lagrange Multiplier
37	_alpha = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );	8✔
38
39	// Initialize kinetic density at grid cells
40	_kineticDensity = VectorVector( _nCells, Vector( _nq, 0.0 ) );	8✔
41
42	// Limiter variables
43	_solDx = VectorVector( _nCells, Vector( _nq, 0.0 ) );	8✔
44	_solDy = VectorVector( _nCells, Vector( _nq, 0.0 ) );	8✔
45	_limiter = VectorVector( _nCells, Vector( _nq, 0.0 ) );	8✔
46
47	// Initialize and Pre-Compute Moments at quadrature points
48	_momentBasis = VectorVector( _nq, Vector( _nSystem, 0.0 ) );	8✔
49	ComputeMoments();	8✔
50
51	// Initialize Scatter Matrix --
52	_scatterMatDiag = Vector( _nSystem, 0.0 );	8✔
53	ComputeScatterMatrix();	8✔
54	}	8✔
55
56	MNSolver::~MNSolver() {	16✔
57	delete _entropy;	8✔
58	delete _optimizer;	8✔
59	delete _basis;	8✔
60	}	16✔
61		8✔
62	void MNSolver::ComputeScatterMatrix() {
63
64	// --- Isotropic ---
65	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	8✔
66	_scatterMatDiag[0] = 1.0;	8✔
67	for( unsigned idx_diag = 1; idx_diag < _nSystem; idx_diag++ ) {	8✔
68	_scatterMatDiag[idx_diag] = 0.0; // SPHERICAL_HARMONICS are orthogonal basis	8✔
69	}	8✔
70	}	8✔
71	else { // SPHERICAL_MONOMIALS
72	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {	8✔
73	_scatterMatDiag[idx_sys] = 0.0;
74	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
75	_scatterMatDiag[idx_sys] += _momentBasis[idx_quad][idx_sys] * _weights[idx_quad];	8✔
76	}	8✔
77	if( _settings->GetDim() == 1 ) {	72✔
78	_scatterMatDiag[idx_sys] /= 2;	64✔
79	}
80	else if( _settings->GetDim() == 2 ) {
81	_scatterMatDiag[idx_sys] /= M_PI;
82	}	×
83	else { // 3D	×
84	_scatterMatDiag[idx_sys] /= 4.0 * M_PI;	×
85	}	×
86	}
87	}	×
88	}	×
89
90	void MNSolver::ComputeMoments() {	×
91	double my, phi;	×
92
93	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
94	my = _quadPointsSphere[idx_quad][0];	×
95	phi = _quadPointsSphere[idx_quad][1];
96	_momentBasis[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
97	}
98	}	8✔
99
100	void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {	8✔
101	// double entropyReconstruction = 0.0;
102
103	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) { // reset solution	1,032✔
104	_sol[idx_cell][idx_sys] = 0.0;	1,024✔
105	}	1,024✔
106	_optimizer->ReconstructMoments( _sol[idx_cell], _alpha[idx_cell], _momentBasis );	1,024✔
107	}
108		8✔
109	void MNSolver::IterPreprocessing( unsigned /idx_pseudotime/ ) {
110		12,288✔
111	// ------- Entropy closure Step ----------------
112	Vector alpha_norm_dummy( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE
113		122,880✔
114	_optimizer->SolveMultiCell( _alpha, _sol, _momentBasis, alpha_norm_dummy ); // parallel	110,592✔
115
116	// ------- Solution reconstruction step ----	12,288✔
117	#pragma omp parallel for	12,288✔
118	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
119	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {	28✔
120	// compute the kinetic density at all grid cells
121	_kineticDensity[idx_cell][idx_quad] = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _momentBasis[idx_quad] ) );
122	}	56✔
123	if( _settings->GetRealizabilityReconstruction() ) ComputeRealizableSolution( idx_cell );
124	}	28✔
125	// ------ Compute slope limiters and cell gradients ---
126	if( _reconsOrder > 1 ) {
127	_mesh->ComputeSlopes( _nq, _solDx, _solDy, _kineticDensity ); // parallel	28✔
128	_mesh->ComputeLimiter( _nq, _solDx, _solDy, _kineticDensity, _limiter ); // parallel
129	}
130	}
131
132	void MNSolver::ComputeScalarFlux() {
133	double firstMomentScaleFactor = 4 * M_PI;
134	if( _settings->GetProblemName() == PROBLEM_Aircavity1D \|\| _settings->GetProblemName() == PROBLEM_Linesource1D \|\|
135	_settings->GetProblemName() == PROBLEM_Checkerboard1D ) {
136	firstMomentScaleFactor = 2.0;	28✔
137	}	6✔
138	if( _settings->GetDim() == 2 ) {	6✔
139	firstMomentScaleFactor = M_PI;
140	}	28✔
141	#pragma omp parallel for
142	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	28✔
143	_scalarFluxNew[idx_cell] = _sol[idx_cell][0] * firstMomentScaleFactor;	28✔
144	}	56✔
145	}	28✔
146		×
147	void MNSolver::FluxUpdate() {
148	// Loop over all spatial cells	28✔
149	if( _settings->GetProblemName() == PROBLEM_Aircavity1D \|\| _settings->GetProblemName() == PROBLEM_Linesource1D \|\|	×
150	_settings->GetProblemName() == PROBLEM_Checkerboard1D ) {
151	FluxUpdatePseudo1D();	28✔
152	}
153	else {
154	FluxUpdatePseudo2D();
155	}	28✔
156	}
157		28✔
158	void MNSolver::FluxUpdatePseudo1D() {
159	// Loop over the grid cells	56✔
160	#pragma omp parallel for	28✔
161	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {	×
162	// Dirichlet Boundaries stayd
163	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
164	double solL, solR, kineticFlux;	28✔
165	Vector flux( _nSystem, 0.0 );
166	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {	28✔
167	kineticFlux = 0.0; // reset temorary flux
168		×
169	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); idx_nbr++ ) {
170	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {	×
171	// Boundary cells are first order and mirror ghost cells
172	solL = _kineticDensity[idx_cell][idx_quad];
173	solR = solL;
174	}
175	else {
176	// interior cell
177	unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
178	if( _reconsOrder == 1 ) {
179	solL = _kineticDensity[idx_cell][idx_quad];
180	solR = _kineticDensity[_neighbors[idx_cell][idx_nbr]][idx_quad];
181	}
182	else if( _reconsOrder == 2 ) {
183	solL = _kineticDensity[idx_cell][idx_quad] +
184	_limiter[idx_cell][idx_quad] *
185	( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
186	_solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
187	solR = _kineticDensity[nbr_glob][idx_quad] +
188	_limiter[nbr_glob][idx_quad] *
189	( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
190	_solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
191	}
192	else {
193	ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION );
194	}
195	}
196	// Kinetic flux
197	kineticFlux += _g->Flux1D( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
198	}
199	// Solution flux
200	flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
201	}
202	_solNew[idx_cell] = flux;
203	}
204	}
205
206	void MNSolver::FluxUpdatePseudo2D() {
207	// Loop over the grid cells
208	#pragma omp parallel for
209	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
210	// Dirichlet Boundaries stayd
211	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
212	double solL, solR, kineticFlux;
213	Vector flux( _nSystem, 0.0 );
214	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
215	kineticFlux = 0.0; // reset temorary flux
216	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); idx_nbr++ ) {	28✔
217	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
218	// Boundary cells are first order and mirror ghost cells	28✔
219	solL = _kineticDensity[idx_cell][idx_quad];
220	solR = solL;
221	}
222	else {
223	// interior cell
224	unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
225	if( _reconsOrder == 1 ) {
226	solL = _kineticDensity[idx_cell][idx_quad];
227	solR = _kineticDensity[_neighbors[idx_cell][idx_nbr]][idx_quad];
228	}
229	else if( _reconsOrder == 2 ) {
230	solL = _kineticDensity[idx_cell][idx_quad] +
231	_limiter[idx_cell][idx_quad] *
232	( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
233	_solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
234	solR = _kineticDensity[nbr_glob][idx_quad] +
235	_limiter[nbr_glob][idx_quad] *
236	( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
237	_solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
238	}
239	else {
240	ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION );
241	}
242	}
243	// Kinetic flux
244	kineticFlux += _g->Flux( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
245	}
246	// Solution flux
247	flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
248	}
249	_solNew[idx_cell] = flux;
250	}
251	}
252
253	void MNSolver::FVMUpdate( unsigned idx_iter ) {
254	if( _Q.size() == 1u && _sigmaT.size() == 1u && _sigmaS.size() == 1u ) { // Physics constant in time
255	// Loop over the grid cells
256	#pragma omp parallel for
257	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
258	// Dirichlet Boundaries stay
259	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
260	// Flux update
261	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {	28✔
262	_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] /* old solution */
263	- ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */	28✔
264	+ _dT * _sigmaS[0][idx_cell] * _sol[idx_cell][0] * _scatterMatDiag[idx_sys] /* scattering gain */	28✔
265	- _dT * ( _sigmaT[0][idx_cell] ) * _sol[idx_cell][idx_sys]; /* scattering and absorbtion loss */
266	}	28✔
267	// Source Term
268	_solNew[idx_cell] += _dT * _Q[0][idx_cell];
269	}
270	}
271	else {
272	// Loop over the grid cells
273	#pragma omp parallel for
274	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
275	// Dirichlet Boundaries stay
276	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
277	// Flux update
278	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
279	_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] /* old solution */
280	- ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
281	+ _dT * _sigmaS[idx_iter][idx_cell] * _sol[idx_cell][0] * _scatterMatDiag[idx_sys] /* scattering gain */
282	- _dT * ( _sigmaT[idx_iter][idx_cell] ) * _sol[idx_cell][idx_sys]; /* scattering and absorbtion loss */
NEW 283	}	×
284	// Source Term
285	_solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell];
286	}
287	}
288	}
289
290	void MNSolver::PrepareVolumeOutput() {
291	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
292	_outputFieldNames.resize( nGroups );
293	_outputFields.resize( nGroups );
294	// Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
295	#pragma omp parallel for
296	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
297	// Prepare all Output Fields per group
298	// Different procedure, depending on the Group...	28✔
299	switch( _settings->GetVolumeOutput()[idx_group] ) {
300	case MINIMAL:	8✔
301	// Currently only one entry ==> rad flux	8✔
302	_outputFields[idx_group].resize( 1 );	8✔
303	_outputFieldNames[idx_group].resize( 1 );	8✔
304	_outputFields[idx_group][0].resize( _nCells );
305	_outputFieldNames[idx_group][0] = "radiation flux density";	8✔
306	break;
307	case MOMENTS:
308	// As many entries as there are moments in the system
309	_outputFields[idx_group].resize( _nSystem );
310	_outputFieldNames[idx_group].resize( _nSystem );
311	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
312	if( _settings->GetDim() == 3 ) {
313	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) { // 3D
314	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
315	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
316	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
317	std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
318	}
319	}
320	}
321	else if( _settings->GetDim() == 2 ) {
322	unsigned count = 0;
323	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
324	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
325	if( ( idx_l + idx_k ) % 2 == 0 ) {
326	_outputFields[idx_group][count].resize( _nCells );
327	_outputFieldNames[idx_group][count] =
328	std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
329	count++;
330	}
331	}
332	}
333	}
334	else if( _settings->GetDim() == 1 ) {
335	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
336	_outputFields[idx_group][idx_l].resize( _nCells );
337	_outputFieldNames[idx_group][idx_l] = std::string( "u_" + std::to_string( idx_l ) + "^0" );
338	}
339	}
340	}
341	else {
342	for( unsigned idx_l = 0; idx_l <= _polyDegreeBasis; idx_l++ ) {
343	unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
344	for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
345	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
346	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
347	std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
348	}
349	}
350	}
351	break;
352	case DUAL_MOMENTS:
353	// As many entries as there are moments in the system
354	_outputFields[idx_group].resize( _nSystem );
355	_outputFieldNames[idx_group].resize( _nSystem );
356	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
357	if( _settings->GetDim() == 3 ) {
358	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
359	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
360	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
361	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
362	std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
363	}
364	}
365	}
366	else if( _settings->GetDim() == 2 ) {
367	unsigned count = 0;
368	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
369	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
370	if( ( idx_l + idx_k ) % 2 == 0 ) {
371	_outputFields[idx_group][count].resize( _nCells );
372	_outputFieldNames[idx_group][count] =
373	std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
374	count++;
375	}
376	}
377	}
378	}
379	else if( _settings->GetDim() == 1 ) {
380	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
381	_outputFields[idx_group][idx_l].resize( _nCells );
382	_outputFieldNames[idx_group][idx_l] = std::string( "alpha_" + std::to_string( idx_l ) + "^0" );
383	}
384	}
385	}
386	else { // SPHERICAL_MONOMIALS
387	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
388	unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
389	for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
390	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
391	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
392	std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
393	}
394	}
395	}
396	break;
397	case ANALYTIC:
398	// one entry per cell
399	_outputFields[idx_group].resize( 1 );
400	_outputFieldNames[idx_group].resize( 1 );
401	_outputFields[idx_group][0].resize( _nCells );
402	_outputFieldNames[idx_group][0] = std::string( "analytic radiation flux density" );
403	break;
404	default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
405	}
406	}
407	}
408
409	void MNSolver::WriteVolumeOutput( unsigned idx_iter ) {
410	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
411	// Check if volume output fields are written to file this iteration
412	if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) \|\|
413	( idx_iter == _nEnergies - 1 ) /* need sol at last iteration */ ) {
414	#pragma omp parallel for
415	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
416	switch( _settings->GetVolumeOutput()[idx_group] ) {
417	case MINIMAL:	8✔
418	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
419	_outputFields[idx_group][0][idx_cell] = _scalarFluxNew[idx_cell];	28✔
420	}	28✔
421	break;
422	case MOMENTS:	56✔
423	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {	28✔
424	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	8✔
425	_outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
426	}
427	}
428	break;
429	case DUAL_MOMENTS:
430	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
431	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
432	_outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];
433	}
434	}
435	break;
436	case ANALYTIC:
437	// Compute total "mass" of the system ==> to check conservation properties
438	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
439	double time = idx_iter * _dT;
440	_outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
441	_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_iter][idx_cell] );
442	}
443	break;
444	default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
445	}
446	}
447	}
448	}

CSMMLab / KiT-RT / #95

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous