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

0.0

/src/solvers/csdmnsolver.cpp

#include "solvers/csdmnsolver.hpp"
#include "common/config.hpp"
#include "common/mesh.hpp"
#include "entropies/entropybase.hpp"
#include "fluxes/numericalflux.hpp"
#include "optimizers/optimizerbase.hpp"
#include "solvers/csdpn_starmap_constants.hpp"
#include "toolboxes/interpolation.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"
#include "velocitybasis/sphericalbase.hpp"

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

CSDMNSolver::CSDMNSolver( Config* settings ) : MNSolver( settings ) {
    // --- Initialize Dose ---
    _dose = std::vector<double>( _nCells, 0.0 );

    // --- Compute transformed energy grid ---

    // determine transformed energy grid for tabulated grid
    Vector E_transformed( E_trans.size(), 0.0 );
    for( unsigned i = 1; i < E_trans.size(); ++i )
        E_transformed[i] = E_transformed[i - 1] + ( E_tab[i] - E_tab[i - 1] ) / 2 * ( 1.0 / S_tab[i] + 1.0 / S_tab[i - 1] );

    // determine minimal and maximal energies
    double minE = 5e-5;
    double maxE = _settings->GetMaxEnergyCSD();

    // define interpolation from energies to corresponding transformed energies \tilde{E} (without substraction of eMaxTrafo)
    Interpolation interpEToTrafo( E_tab, E_transformed );
    double eMaxTrafo = interpEToTrafo( maxE );
    double eMinTrafo = interpEToTrafo( minE );

    // check what happens if we intepolate back
    Interpolation interpTrafoToE( E_transformed, E_tab );

    // define linear grid in fully transformed energy \tilde\tilde E (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
    _eTrafo = blaze::linspace( _nEnergies, eMaxTrafo - eMaxTrafo, eMaxTrafo - eMinTrafo );
    // double dETrafo = _eTrafo[1] - _eTrafo[0];

    // compute Corresponding original energies
    for( unsigned n = 0; n < _nEnergies; ++n ) {
        _energies[n] = interpTrafoToE( eMaxTrafo - _eTrafo[n] );
    }

    // --- evaluate corresponding stopping powers and transport coefficients

    // compute stopping powers
    Vector etmp = E_tab;
    Vector stmp = S_tab;
    Interpolation interpS( etmp, stmp );
    _sigmaTAtEnergy = Vector( _polyDegreeBasis + 1, 0.0 );

    // compute stopping power between energies for dose computation
    double dE = _eTrafo[2] - _eTrafo[1];
    Vector eTrafoMid( _nEnergies - 1 );
    for( unsigned n = 0; n < _nEnergies - 1; ++n ) {
        eTrafoMid[n] = _eTrafo[n] + dE / 2;
    }
    // compute Corresponding original energies at intermediate points
    Vector eMid( _nEnergies - 1 );
    for( unsigned n = 0; n < _nEnergies - 1; ++n ) {
        eMid[n] = interpTrafoToE( eMaxTrafo - eTrafoMid[n] );
    }
    _sMid = interpS( eMid );
}

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, alpha_norm_dummy );

    // ------- 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 density normalized slope limiters and cell gradients ---
    if( _reconsOrder > 1 ) {
        VectorVector solDivRho = _sol;
        for( unsigned j = 0; j < _nCells; ++j ) {
            solDivRho[j] = _kineticDensity[j] / _density[j];
        }
        _mesh->ComputeSlopes( _nq, _solDx, _solDy, solDivRho );
        _mesh->ComputeLimiter( _nq, _solDx, _solDy, solDivRho, _limiter );
    }

    // ------ evaluate scatter coefficient at current energy level
    Vector sigmaSAtEnergy( _polyDegreeBasis + 1, 0.0 );
    // compute scattering cross section at current energy
    for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {
        // setup interpolation from E to sigma at degree idx_degree
        Interpolation interp( E_ref, blaze::column( sigma_ref, idx_degree ) );
        sigmaSAtEnergy[idx_degree] = interp( _energies[idx_iter + 1] );
    }
    for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {
        _sigmaTAtEnergy[idx_degree] = ( sigmaSAtEnergy[0] - sigmaSAtEnergy[idx_degree] );
    }
}

void CSDMNSolver::SolverPreprocessing() {
    // cross sections do not need to be transformed to ETilde energy grid since e.g. TildeSigmaT(ETilde) = SigmaT(E(ETilde))
}

void CSDMNSolver::IterPostprocessing( unsigned idx_iter ) {
    // --- Compute Flux for solution and Screen Output ---
    ComputeScalarFlux();

    // --- Compute Dose ---
#pragma omp parallel for
    for( unsigned j = 0; j < _nCells; ++j ) {
        if( idx_iter > 0 && idx_iter < _nEnergies - 1 ) {
            _dose[j] += _dT * ( _sol[j][0] * _sMid[idx_iter] ) / _density[j];    // update dose with trapezoidal rule // diss Kerstin
        }
        else {
            _dose[j] += 0.5 * _dT * ( _sol[j][0] * _sMid[idx_iter] ) / _density[j];
        }
    }
}

void CSDMNSolver::FluxUpdate() {
// 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;

        //--- Integration of moments of flux ---
        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] / _density[idx_cell];
                    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] / _density[idx_cell];
                        solR = _kineticDensity[nbr_glob][idx_quad] / _density[nbr_glob];
                    }
                    else if( _reconsOrder == 2 ) {
                        solL = _kineticDensity[idx_cell][idx_quad] / _density[idx_cell] +
                               _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] / _density[nbr_glob] +
                               _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->FluxXZ( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
            }
            // Solution flux
            flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
        }

        _solNew[idx_cell] = flux;    // ConstructFlux( idx_cell );
    }
}

void CSDMNSolver::FVMUpdate( unsigned idx_energy ) {
    bool implicitScattering = true;
    // transform energy difference
    _dT = fabs( _eTrafo[idx_energy + 1] - _eTrafo[idx_energy] );
    // loop over all spatial cells
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        //  Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                int idx_sys = _basis->GetGlobalIndexBasis( idx_l, idx_k );
                if( implicitScattering ) {
                    _solNew[idx_cell][idx_sys] =
                        _sol[idx_cell][idx_sys] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys];            /* cell averaged flux */
                    _solNew[idx_cell][idx_sys] = _solNew[idx_cell][idx_sys] / ( 1.0 + _dT * _sigmaTAtEnergy[idx_l] ); /* implicit scattering */
                }
                else {
                    _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] -
                                                 ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys]   /* cell averaged flux */
                                                 - _dT * _sol[idx_cell][idx_sys] * _sigmaTAtEnergy[idx_l]; /* scattering */
                }
            }
            // Source Term TODO
            // _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
        }
    }
}

void CSDMNSolver::PrepareVolumeOutput() {
    // std::cout << "Prepare Volume Output...";
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();

    _outputFieldNames.resize( nGroups );
    _outputFields.resize( nGroups );

    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
        // Prepare all Output Fields per group
        unsigned glob_idx = 0;

        // 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 MEDICAL:
                _outputFields[idx_group].resize( 2 );
                _outputFieldNames[idx_group].resize( 2 );

                // Dose
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "dose";
                // Normalized Dose
                _outputFields[idx_group][1].resize( _nCells );
                _outputFieldNames[idx_group][1] = "normalized dose";
                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 ) {
                    for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
                        for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                            glob_idx = _basis->GetGlobalIndexBasis( idx_l, idx_k );
                            _outputFields[idx_group][glob_idx].resize( _nCells );
                            _outputFieldNames[idx_group][glob_idx] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
                    }
                }
                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 ) {
                    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 {
                    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;
            default: ErrorMessages::Error( "Volume Output Group not defined for CSD MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

void CSDMNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
    double maxDose;
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
        ( idx_pseudoTime == _nIter - 1 ) /* need sol at last iteration */ ) {

        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 MEDICAL:
                    // Compute Dose
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        _outputFields[idx_group][0][idx_cell] = _dose[idx_cell];
                    }
                    // Compute normalized dose
                    _outputFields[idx_group][1] = _outputFields[idx_group][0];

                    maxDose = *std::max_element( _outputFields[idx_group][0].begin(), _outputFields[idx_group][0].end() );

                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        _outputFields[idx_group][1][idx_cell] /= maxDose;
                    }
                    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;
                default: ErrorMessages::Error( "Volume Output Group not defined for CSD MN Solver!", CURRENT_FUNCTION ); break;
            }
        }
    }
    if( idx_pseudoTime == _nEnergies - 2 ) {
        std::ofstream out( _settings->GetOutputFile().append( ".txt" ) );
        unsigned nx = _settings->GetNCells();

        for( unsigned j = 0; j < nx; ++j ) {
            out << _cellMidPoints[j][0] << " " << _cellMidPoints[j][1] << " " << _dose[j] << std::endl;
        }
        out.close();
    }
}

double CSDMNSolver::NormPDF( double x, double mu, double sigma ) {
    return INV_SQRT_2PI / sigma * std::exp( -( ( x - mu ) * ( x - mu ) ) / ( 2.0 * sigma * sigma ) );
}

1	#include "solvers/csdmnsolver.hpp"
2	#include "common/config.hpp"
3	#include "common/mesh.hpp"
4	#include "entropies/entropybase.hpp"
5	#include "fluxes/numericalflux.hpp"
6	#include "optimizers/optimizerbase.hpp"
7	#include "solvers/csdpn_starmap_constants.hpp"
8	#include "toolboxes/interpolation.hpp"
9	#include "toolboxes/textprocessingtoolbox.hpp"
10	#include "velocitybasis/sphericalbase.hpp"
11
12	// externals
13	#include "spdlog/spdlog.h"
14
15	CSDMNSolver::CSDMNSolver( Config* settings ) : MNSolver( settings ) {	×
16	// --- Initialize Dose ---
17	_dose = std::vector<double>( _nCells, 0.0 );	×
18
19	// --- Compute transformed energy grid ---
20
21	// determine transformed energy grid for tabulated grid
22	Vector E_transformed( E_trans.size(), 0.0 );	×
23	for( unsigned i = 1; i < E_trans.size(); ++i )	×
24	E_transformed[i] = E_transformed[i - 1] + ( E_tab[i] - E_tab[i - 1] ) / 2 * ( 1.0 / S_tab[i] + 1.0 / S_tab[i - 1] );	×
25
26	// determine minimal and maximal energies
27	double minE = 5e-5;	×
28	double maxE = _settings->GetMaxEnergyCSD();	×
29
30	// define interpolation from energies to corresponding transformed energies \tilde{E} (without substraction of eMaxTrafo)
31	Interpolation interpEToTrafo( E_tab, E_transformed );	×
32	double eMaxTrafo = interpEToTrafo( maxE );	×
33	double eMinTrafo = interpEToTrafo( minE );	×
34
35	// check what happens if we intepolate back
36	Interpolation interpTrafoToE( E_transformed, E_tab );	×
37
38	// define linear grid in fully transformed energy \tilde\tilde E (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
39	_eTrafo = blaze::linspace( _nEnergies, eMaxTrafo - eMaxTrafo, eMaxTrafo - eMinTrafo );	×
40	// double dETrafo = _eTrafo[1] - _eTrafo[0];
41
42	// compute Corresponding original energies
43	for( unsigned n = 0; n < _nEnergies; ++n ) {	×
44	_energies[n] = interpTrafoToE( eMaxTrafo - _eTrafo[n] );	×
45	}
46
47	// --- evaluate corresponding stopping powers and transport coefficients
48
49	// compute stopping powers
50	Vector etmp = E_tab;	×
51	Vector stmp = S_tab;	×
52	Interpolation interpS( etmp, stmp );	×
53	_sigmaTAtEnergy = Vector( _polyDegreeBasis + 1, 0.0 );	×
54
55	// compute stopping power between energies for dose computation
56	double dE = _eTrafo[2] - _eTrafo[1];	×
57	Vector eTrafoMid( _nEnergies - 1 );	×
58	for( unsigned n = 0; n < _nEnergies - 1; ++n ) {	×
59	eTrafoMid[n] = _eTrafo[n] + dE / 2;	×
60	}
61	// compute Corresponding original energies at intermediate points
62	Vector eMid( _nEnergies - 1 );	×
63	for( unsigned n = 0; n < _nEnergies - 1; ++n ) {	×
64	eMid[n] = interpTrafoToE( eMaxTrafo - eTrafoMid[n] );	×
65	}
66	_sMid = interpS( eMid );	×
67	}
68
69	CSDMNSolver::~CSDMNSolver() {}	×
70		×
71	void CSDMNSolver::IterPreprocessing( unsigned idx_iter ) {	×
72
73	// ------- Entropy closure Step ----------------	×
74	Vector alpha_norm_dummy( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE
75
76	_optimizer->SolveMultiCell( _alpha, _sol, _momentBasis, alpha_norm_dummy );	×
77
78	// ------- Solution reconstruction step ----	×
79	#pragma omp parallel for
80	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
81	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {	×
82	// compute the kinetic density at all grid cells
83	_kineticDensity[idx_cell][idx_quad] = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _momentBasis[idx_quad] ) );
84	}
85	if( _settings->GetRealizabilityReconstruction() ) ComputeRealizableSolution( idx_cell );
86	}
87
88	// ------ Compute density normalized slope limiters and cell gradients ---
89	if( _reconsOrder > 1 ) {
90	VectorVector solDivRho = _sol;
91	for( unsigned j = 0; j < _nCells; ++j ) {	×
92	solDivRho[j] = _kineticDensity[j] / _density[j];	×
93	}	×
94	_mesh->ComputeSlopes( _nq, _solDx, _solDy, solDivRho );	×
95	_mesh->ComputeLimiter( _nq, _solDx, _solDy, solDivRho, _limiter );
96	}	×
97		×
98	// ------ evaluate scatter coefficient at current energy level
99	Vector sigmaSAtEnergy( _polyDegreeBasis + 1, 0.0 );
100	// compute scattering cross section at current energy
101	for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {	×
102	// setup interpolation from E to sigma at degree idx_degree
103	Interpolation interp( E_ref, blaze::column( sigma_ref, idx_degree ) );	×
104	sigmaSAtEnergy[idx_degree] = interp( _energies[idx_iter + 1] );
105	}	×
106	for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {	×
107	_sigmaTAtEnergy[idx_degree] = ( sigmaSAtEnergy[0] - sigmaSAtEnergy[idx_degree] );
108	}	×
109	}	×
110
111	void CSDMNSolver::SolverPreprocessing() {
112	// cross sections do not need to be transformed to ETilde energy grid since e.g. TildeSigmaT(ETilde) = SigmaT(E(ETilde))
113	}	×
114
115	void CSDMNSolver::IterPostprocessing( unsigned idx_iter ) {
116	// --- Compute Flux for solution and Screen Output ---
NEW 117	ComputeScalarFlux();	×
118
119	// --- Compute Dose ---	×
120	#pragma omp parallel for
121	for( unsigned j = 0; j < _nCells; ++j ) {
122	if( idx_iter > 0 && idx_iter < _nEnergies - 1 ) {	×
123	_dose[j] += _dT * ( _sol[j][0] * _sMid[idx_iter] ) / _density[j]; // update dose with trapezoidal rule // diss Kerstin
124	}
125	else {
126	_dose[j] += 0.5 * _dT * ( _sol[j][0] * _sMid[idx_iter] ) / _density[j];
127	}
128	}
129	}
130
131	void CSDMNSolver::FluxUpdate() {
132	// Loop over the grid cells
133	#pragma omp parallel for	×
134	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
135	// Dirichlet Boundaries stayd	×
136	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
137
138	//--- Integration of moments of flux ---
139	double solL, solR, kineticFlux;
140	Vector flux( _nSystem, 0.0 );
141
142	for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
143	kineticFlux = 0.0; // reset temorary flux
144
145	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); idx_nbr++ ) {
146	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
147	// Boundary cells are first order and mirror ghost cells
148	solL = _kineticDensity[idx_cell][idx_quad] / _density[idx_cell];
149	solR = solL;
150	}
151	else {
152	// interior cell
153	unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
154	if( _reconsOrder == 1 ) {
155	solL = _kineticDensity[idx_cell][idx_quad] / _density[idx_cell];
156	solR = _kineticDensity[nbr_glob][idx_quad] / _density[nbr_glob];
157	}
158	else if( _reconsOrder == 2 ) {
159	solL = _kineticDensity[idx_cell][idx_quad] / _density[idx_cell] +
160	_limiter[idx_cell][idx_quad] *
161	( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
162	_solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
163	solR = _kineticDensity[nbr_glob][idx_quad] / _density[nbr_glob] +
164	_limiter[nbr_glob][idx_quad] *
165	( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
166	_solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
167	}
168	else {
169	ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION );
170	}
171	}
172	// Kinetic flux
173	kineticFlux += _g->FluxXZ( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
174	}
175	// Solution flux
176	flux += _momentBasis[idx_quad] * ( _weights[idx_quad] * kineticFlux );
177	}
178
179	_solNew[idx_cell] = flux; // ConstructFlux( idx_cell );
180	}
181	}
182
183	void CSDMNSolver::FVMUpdate( unsigned idx_energy ) {
184	bool implicitScattering = true;
185	// transform energy difference	×
NEW 186	_dT = fabs( _eTrafo[idx_energy + 1] - _eTrafo[idx_energy] );	×
187	// loop over all spatial cells
188	#pragma omp parallel for	×
189	for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
190	// Dirichlet cells stay at IC, farfield assumption	×
191	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
192	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
193	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
194	int idx_sys = _basis->GetGlobalIndexBasis( idx_l, idx_k );
195	if( implicitScattering ) {
196	_solNew[idx_cell][idx_sys] =
197	_sol[idx_cell][idx_sys] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys]; /* cell averaged flux */
198	_solNew[idx_cell][idx_sys] = _solNew[idx_cell][idx_sys] / ( 1.0 + _dT * _sigmaTAtEnergy[idx_l] ); /* implicit scattering */
199	}
200	else {
201	_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] -
202	( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
203	- _dT * _sol[idx_cell][idx_sys] * _sigmaTAtEnergy[idx_l]; /* scattering */
204	}
205	}
206	// Source Term TODO
207	// _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
208	}
209	}
210	}
211
212	void CSDMNSolver::PrepareVolumeOutput() {
213	// std::cout << "Prepare Volume Output...";
214	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();	×
215
216	_outputFieldNames.resize( nGroups );	×
217	_outputFields.resize( nGroups );
218		×
219	// Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT	×
220	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
221	// Prepare all Output Fields per group
222	unsigned glob_idx = 0;	×
223
224	// Different procedure, depending on the Group...	×
225	switch( _settings->GetVolumeOutput()[idx_group] ) {
226	case MINIMAL:
227	// Currently only one entry ==> rad flux	×
228	_outputFields[idx_group].resize( 1 );	×
229	_outputFieldNames[idx_group].resize( 1 );
230		×
231	_outputFields[idx_group][0].resize( _nCells );	×
232	_outputFieldNames[idx_group][0] = "radiation flux density";
233	break;	×
234		×
235	case MEDICAL:	×
236	_outputFields[idx_group].resize( 2 );
237	_outputFieldNames[idx_group].resize( 2 );	×
238		×
239	// Dose	×
240	_outputFields[idx_group][0].resize( _nCells );
241	_outputFieldNames[idx_group][0] = "dose";
242	// Normalized Dose	×
243	_outputFields[idx_group][1].resize( _nCells );	×
244	_outputFieldNames[idx_group][1] = "normalized dose";
245	break;	×
246	case MOMENTS:	×
247	// As many entries as there are moments in the system	×
248	_outputFields[idx_group].resize( _nSystem );	×
249	_outputFieldNames[idx_group].resize( _nSystem );
250	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
251	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {	×
252	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {	×
253	glob_idx = _basis->GetGlobalIndexBasis( idx_l, idx_k );	×
254	_outputFields[idx_group][glob_idx].resize( _nCells );	×
255	_outputFieldNames[idx_group][glob_idx] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );	×
256	}	×
257	}	×
258	}
259	else {
260	for( unsigned idx_l = 0; idx_l <= _polyDegreeBasis; idx_l++ ) {
261	unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
262	for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {	×
263	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );	×
264	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =	×
265	std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );	×
266	}	×
267	}	×
268	}
269	break;
270	case DUAL_MOMENTS:
271	// As many entries as there are moments in the system	×
272	_outputFields[idx_group].resize( _nSystem );	×
273	_outputFieldNames[idx_group].resize( _nSystem );
274		×
275	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
276	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
277	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {	×
278	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );	×
279	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =	×
280	std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );	×
281	}	×
282	}	×
283	}
284	else {
285	for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
286	unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
287	for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {	×
288	_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );	×
289	_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =	×
290	std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );	×
291	}	×
292	}	×
293	}
294	break;
295	default: ErrorMessages::Error( "Volume Output Group not defined for CSD MN Solver!", CURRENT_FUNCTION ); break;
296	}	×
297	}	×
298	}
299
300	void CSDMNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
301	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
302	double maxDose;	×
303	if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) \|\|	×
304	( idx_pseudoTime == _nIter - 1 ) /* need sol at last iteration */ ) {
305		×
306	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {	×
307	switch( _settings->GetVolumeOutput()[idx_group] ) {
308	case MINIMAL:	×
309	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
NEW 310	_outputFields[idx_group][0][idx_cell] = _scalarFluxNew[idx_cell];	×
311	}	×
312	break;	×
313
314	case MEDICAL:	×
315	// Compute Dose
316	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
317	_outputFields[idx_group][0][idx_cell] = _dose[idx_cell];
318	}	×
319	// Compute normalized dose	×
320	_outputFields[idx_group][1] = _outputFields[idx_group][0];
321
322	maxDose = *std::max_element( _outputFields[idx_group][0].begin(), _outputFields[idx_group][0].end() );	×
323
324	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
325	_outputFields[idx_group][1][idx_cell] /= maxDose;
326	}	×
327	break;	×
328	case MOMENTS:
329	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {	×
330	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
331	_outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];	×
332	}	×
333	}	×
334	break;
335	case DUAL_MOMENTS:
336	for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {	×
337	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
338	_outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];	×
339	}	×
340	}	×
341	break;
342	default: ErrorMessages::Error( "Volume Output Group not defined for CSD MN Solver!", CURRENT_FUNCTION ); break;
343	}	×
344	}	×
345	}
346	if( idx_pseudoTime == _nEnergies - 2 ) {
347	std::ofstream out( _settings->GetOutputFile().append( ".txt" ) );
348	unsigned nx = _settings->GetNCells();	×
349		×
350	for( unsigned j = 0; j < nx; ++j ) {	×
351	out << _cellMidPoints[j][0] << " " << _cellMidPoints[j][1] << " " << _dose[j] << std::endl;
352	}	×
353	out.close();	×
354	}
355	}	×
356
357	double CSDMNSolver::NormPDF( double x, double mu, double sigma ) {
358	return INV_SQRT_2PI / sigma * std::exp( -( ( x - mu ) * ( x - mu ) ) / ( 2.0 * sigma * sigma ) );
359	}	×

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