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/csdsnsolver.cpp

#include "solvers/csdsnsolver.hpp"
#include "common/config.hpp"
#include "common/io.hpp"
#include "common/mesh.hpp"
#include "fluxes/numericalflux.hpp"
#include "kernels/scatteringkernelbase.hpp"
#include "problems/problembase.hpp"
#include "quadratures/qproduct.hpp"
#include "quadratures/quadraturebase.hpp"
#include "solvers/csdpn_starmap_constants.hpp"
#include "solvers/solverbase.hpp"
#include "toolboxes/icru.hpp"

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

CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );

    // 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 sigmaT is now done during computation in IterPreprocessing()

    // compute stopping powers
    Vector etmp = E_tab;
    Vector stmp = S_tab;
    Interpolation interpS( etmp, stmp );
    _s = interpS( _energies );

    // 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 );

    _quadPointsSphere  = _quadrature->GetPointsSphere();    // (my,phi,r) gridpoints of the quadrature in spherical cordinates
    SphericalBase* sph = SphericalBase::Create( _settings );
    Matrix Y           = blaze::zero<double>( sph->GetBasisSize(), _nq );

    for( unsigned q = 0; q < _nq; q++ ) {
        blaze::column( Y, q ) = sph->ComputeSphericalBasis( _quadPointsSphere[q][0], _quadPointsSphere[q][1] );
    }
    delete sph;

    _O = blaze::trans( Y );
    _M = _O;    // transposed to simplify weight multiplication. We will transpose _M later again

    unsigned nSph = _M.columns();
    for( unsigned j = 0; j < nSph; ++j ) {
        blaze::column( _M, j ) *= _weights;
    }

    _M.transpose();
    _polyDegreeBasis = settings->GetMaxMomentDegree();
}

// Solver
void CSDSNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
    if( _reconsOrder > 1 ) {
        VectorVector solDivRho = _sol;
        for( unsigned j = 0; j < _nCells; ++j ) {
            solDivRho[j] = _sol[j] / _density[j];
        }
        _mesh->ComputeSlopes( _nq, _solDx, _solDy, solDivRho );
        _mesh->ComputeLimiter( _nq, _solDx, _solDy, solDivRho, _limiter );
    }

    _dT = fabs( _eTrafo[idx_pseudotime + 1] - _eTrafo[idx_pseudotime] );

    Vector sigmaSAtEnergy( _polyDegreeBasis + 1, 0.0 );
    Vector sigmaTAtEnergy( _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_pseudotime + 1] );
    }

    for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {
        sigmaTAtEnergy[idx_degree] = ( sigmaSAtEnergy[0] - sigmaSAtEnergy[idx_degree] );
    }

    // --- Implicit Scattering (First Scatter then Stream - Splitting) ---
#pragma omp parallel for
    for( unsigned j = 0; j < _nCells; ++j ) {
        if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
        Vector u = _M * _sol[j];

        unsigned counter = 0;
        for( int l = 0; l <= (int)_polyDegreeBasis; ++l ) {
            for( int m = -l; m <= l; ++m ) {
                u[counter] = u[counter] / ( 1.0 + _dT * sigmaTAtEnergy[l] );
                counter++;
            }
        }
        _sol[j] = _O * u;
    }
}

void CSDSNSolver::SolverPreprocessing() {}

void CSDSNSolver::IterPostprocessing( unsigned idx_pseudotime ) {
    unsigned n = idx_pseudotime;

    // --- Compute Flux for solution and Screen Output ---
    // ComputeRadFlux(); // do not scale radflux here

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

void CSDSNSolver::FluxUpdate() {
// loop over all spatial cells
#pragma omp parallel for
    for( unsigned j = 0; j < _nCells; ++j ) {
        double solL;
        double solR;
        if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // loop over all ordinates
        for( unsigned i = 0; i < _nq; ++i ) {
            _solNew[j][i] = 0.0;
            // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
            for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
                // store flux contribution on psiNew_sigmaS to save memory
                if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
                    continue;    // adiabatic wall, add nothing
                else {
                    unsigned int nbr_glob = _neighbors[j][idx_neighbor];    // global idx of neighbor cell

                    switch( _reconsOrder ) {
                        case 1: _solNew[j][i] += _g->FluxXZ( _quadPoints[i], _sol[j][i], _sol[nbr_glob][i], _normals[j][idx_neighbor] ); break;
                        // second order solver
                        case 2:
                            // left status of interface
                            solL = _sol[j][i] / _density[j] +
                                   _limiter[j][i] * ( _solDx[j][i] * ( _interfaceMidPoints[j][idx_neighbor][0] - _cellMidPoints[j][0] ) +
                                                      _solDy[j][i] * ( _interfaceMidPoints[j][idx_neighbor][1] - _cellMidPoints[j][1] ) );
                            solR = _sol[nbr_glob][i] / _density[nbr_glob] +
                                   _limiter[nbr_glob][i] *
                                       ( _solDx[nbr_glob][i] * ( _interfaceMidPoints[j][idx_neighbor][0] - _cellMidPoints[nbr_glob][0] ) +
                                         _solDy[nbr_glob][i] * ( _interfaceMidPoints[j][idx_neighbor][1] - _cellMidPoints[nbr_glob][1] ) );

                            // flux evaluation
                            _solNew[j][i] += _g->FluxXZ( _quadPoints[i], solL, solR, _normals[j][idx_neighbor] ) / _areas[j];
                            break;
                            // higher order solver
                        default: ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION ); break;
                    }
                }
            }
        }
    }
}

void CSDSNSolver::FVMUpdate( unsigned /*idx_energy*/ ) {
// loop over all spatial cells
#pragma omp parallel for
    for( unsigned j = 0; j < _nCells; ++j ) {
        if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // loop over all ordinates
        for( unsigned i = 0; i < _nq; ++i ) {
            // time update angular flux with numerical flux and total scattering cross section
            _solNew[j][i] = _sol[j][i] - _dT * _solNew[j][i];
        }
    }
}

// IO
void CSDSNSolver::PrepareVolumeOutput() {
    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

        // 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;

            default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

void CSDSNSolver::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;

                default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO 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();
    }
}

1	#include "solvers/csdsnsolver.hpp"
2	#include "common/config.hpp"
3	#include "common/io.hpp"
4	#include "common/mesh.hpp"
5	#include "fluxes/numericalflux.hpp"
6	#include "kernels/scatteringkernelbase.hpp"
7	#include "problems/problembase.hpp"
8	#include "quadratures/qproduct.hpp"
9	#include "quadratures/quadraturebase.hpp"
10	#include "solvers/csdpn_starmap_constants.hpp"
11	#include "solvers/solverbase.hpp"
12	#include "toolboxes/icru.hpp"
13
14	// externals
15	#include "spdlog/spdlog.h"
16
UNCOV 17	CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {	×
18	_dose = std::vector<double>( _settings->GetNCells(), 0.0 );	×
19
20	// determine transformed energy grid for tabulated grid
21	Vector E_transformed( E_trans.size(), 0.0 );	×
22	for( unsigned i = 1; i < E_trans.size(); ++i )	×
23	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] );	×
24
25	// determine minimal and maximal energies
26	double minE = 5e-5;	×
27	double maxE = _settings->GetMaxEnergyCSD();	×
28
29	// define interpolation from energies to corresponding transformed energies \tilde{E} (without substraction of eMaxTrafo)
30	Interpolation interpEToTrafo( E_tab, E_transformed );	×
31	double eMaxTrafo = interpEToTrafo( maxE );	×
32	double eMinTrafo = interpEToTrafo( minE );	×
33
34	// check what happens if we intepolate back
35	Interpolation interpTrafoToE( E_transformed, E_tab );	×
36
37	// define linear grid in fully transformed energy \tilde\tilde E (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
38	_eTrafo = blaze::linspace( _nEnergies, eMaxTrafo - eMaxTrafo, eMaxTrafo - eMinTrafo );	×
39	// double dETrafo = _eTrafo[1] - _eTrafo[0];
40
41	// compute Corresponding original energies
42	for( unsigned n = 0; n < _nEnergies; ++n ) {	×
43	_energies[n] = interpTrafoToE( eMaxTrafo - _eTrafo[n] );	×
44	}
45
46	// evaluate corresponding stopping powers and transport coefficients
47	// compute sigmaT is now done during computation in IterPreprocessing()
48
49	// compute stopping powers
50	Vector etmp = E_tab;	×
51	Vector stmp = S_tab;	×
52	Interpolation interpS( etmp, stmp );	×
53	_s = interpS( _energies );	×
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	_quadPointsSphere = _quadrature->GetPointsSphere(); // (my,phi,r) gridpoints of the quadrature in spherical cordinates	×
69	SphericalBase* sph = SphericalBase::Create( _settings );	×
70	Matrix Y = blaze::zero<double>( sph->GetBasisSize(), _nq );	×
71
72	for( unsigned q = 0; q < _nq; q++ ) {	×
73	blaze::column( Y, q ) = sph->ComputeSphericalBasis( _quadPointsSphere[q][0], _quadPointsSphere[q][1] );	×
74	}
75	delete sph;	×
76
77	_O = blaze::trans( Y );	×
78	_M = _O; // transposed to simplify weight multiplication. We will transpose _M later again	×
79
80	unsigned nSph = _M.columns();	×
81	for( unsigned j = 0; j < nSph; ++j ) {	×
82	blaze::column( _M, j ) *= _weights;	×
83	}
84
85	_M.transpose();	×
86	_polyDegreeBasis = settings->GetMaxMomentDegree();	×
87	}
88
89	// Solver
90	void CSDSNSolver::IterPreprocessing( unsigned idx_pseudotime ) {	×
91	if( _reconsOrder > 1 ) {	×
92	VectorVector solDivRho = _sol;	×
93	for( unsigned j = 0; j < _nCells; ++j ) {	×
94	solDivRho[j] = _sol[j] / _density[j];	×
95	}
96	_mesh->ComputeSlopes( _nq, _solDx, _solDy, solDivRho );	×
97	_mesh->ComputeLimiter( _nq, _solDx, _solDy, solDivRho, _limiter );	×
98	}
99
NEW 100	_dT = fabs( _eTrafo[idx_pseudotime + 1] - _eTrafo[idx_pseudotime] );	×
101
102	Vector sigmaSAtEnergy( _polyDegreeBasis + 1, 0.0 );	×
103	Vector sigmaTAtEnergy( _polyDegreeBasis + 1, 0.0 );	×
104	// compute scattering cross section at current energy
105	for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {	×
106	// setup interpolation from E to sigma at degree idx_degree
107	Interpolation interp( E_ref, blaze::column( sigma_ref, idx_degree ) );	×
108	sigmaSAtEnergy[idx_degree] = interp( _energies[idx_pseudotime + 1] );	×
109	}
110
111	for( unsigned idx_degree = 0; idx_degree <= _polyDegreeBasis; ++idx_degree ) {	×
112	sigmaTAtEnergy[idx_degree] = ( sigmaSAtEnergy[0] - sigmaSAtEnergy[idx_degree] );	×
113	}
114
115	// --- Implicit Scattering (First Scatter then Stream - Splitting) ---
116	#pragma omp parallel for	×
117	for( unsigned j = 0; j < _nCells; ++j ) {
118	if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
119	Vector u = _M * _sol[j];
120
121	unsigned counter = 0;
122	for( int l = 0; l <= (int)_polyDegreeBasis; ++l ) {
123	for( int m = -l; m <= l; ++m ) {
124	u[counter] = u[counter] / ( 1.0 + _dT * sigmaTAtEnergy[l] );
125	counter++;
126	}
127	}
128	_sol[j] = _O * u;
129	}
130	}
131
NEW 132	void CSDSNSolver::SolverPreprocessing() {}	×
133
134	void CSDSNSolver::IterPostprocessing( unsigned idx_pseudotime ) {	×
135	unsigned n = idx_pseudotime;	×
136
137	// --- Compute Flux for solution and Screen Output ---
138	// ComputeRadFlux(); // do not scale radflux here
139
140	// -- Compute Dose
141	#pragma omp parallel for	×
142	for( unsigned j = 0; j < _nCells; ++j ) {
143	_scalarFluxNew[j] = dot( _sol[j], _weights ); // unscaled rad flux
144	if( n > 0 && n < _nEnergies - 1 ) {
145	_dose[j] += _dT * ( _scalarFluxNew[j] * _sMid[n] ) / _density[j]; // update dose with trapezoidal rule // diss Kerstin
146	}
147	else {
148	_dose[j] += 0.5 * _dT * ( _scalarFluxNew[j] * _sMid[n] ) / _density[j];
149	}
150	}
151	}
152
153	void CSDSNSolver::FluxUpdate() {	×
154	// loop over all spatial cells
155	#pragma omp parallel for	×
156	for( unsigned j = 0; j < _nCells; ++j ) {
157	double solL;
158	double solR;
159	if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
160	// loop over all ordinates
161	for( unsigned i = 0; i < _nq; ++i ) {
162	_solNew[j][i] = 0.0;
163	// loop over all neighbor cells (edges) of cell j and compute numerical fluxes
164	for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
165	// store flux contribution on psiNew_sigmaS to save memory
166	if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
167	continue; // adiabatic wall, add nothing
168	else {
169	unsigned int nbr_glob = _neighbors[j][idx_neighbor]; // global idx of neighbor cell
170
171	switch( _reconsOrder ) {
172	case 1: _solNew[j][i] += _g->FluxXZ( _quadPoints[i], _sol[j][i], _sol[nbr_glob][i], _normals[j][idx_neighbor] ); break;
173	// second order solver
174	case 2:
175	// left status of interface
176	solL = _sol[j][i] / _density[j] +
177	_limiter[j][i] * ( _solDx[j][i] * ( _interfaceMidPoints[j][idx_neighbor][0] - _cellMidPoints[j][0] ) +
178	_solDy[j][i] * ( _interfaceMidPoints[j][idx_neighbor][1] - _cellMidPoints[j][1] ) );
179	solR = _sol[nbr_glob][i] / _density[nbr_glob] +
180	_limiter[nbr_glob][i] *
181	( _solDx[nbr_glob][i] * ( _interfaceMidPoints[j][idx_neighbor][0] - _cellMidPoints[nbr_glob][0] ) +
182	_solDy[nbr_glob][i] * ( _interfaceMidPoints[j][idx_neighbor][1] - _cellMidPoints[nbr_glob][1] ) );
183
184	// flux evaluation
185	_solNew[j][i] += _g->FluxXZ( _quadPoints[i], solL, solR, _normals[j][idx_neighbor] ) / _areas[j];
186	break;
187	// higher order solver
188	default: ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION ); break;
189	}
190	}
191	}
192	}
193	}
194	}
195
196	void CSDSNSolver::FVMUpdate( unsigned /idx_energy/ ) {	×
197	// loop over all spatial cells
198	#pragma omp parallel for	×
199	for( unsigned j = 0; j < _nCells; ++j ) {
200	if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
201	// loop over all ordinates
202	for( unsigned i = 0; i < _nq; ++i ) {
203	// time update angular flux with numerical flux and total scattering cross section
204	_solNew[j][i] = _sol[j][i] - _dT * _solNew[j][i];
205	}
206	}
207	}
208
209	// IO
210	void CSDSNSolver::PrepareVolumeOutput() {	×
211	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();	×
212
213	_outputFieldNames.resize( nGroups );	×
214	_outputFields.resize( nGroups );	×
215
216	// Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
217	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {	×
218	// Prepare all Output Fields per group
219
220	// Different procedure, depending on the Group...
221	switch( _settings->GetVolumeOutput()[idx_group] ) {	×
222	case MINIMAL:	×
223	// Currently only one entry ==> rad flux
224	_outputFields[idx_group].resize( 1 );	×
225	_outputFieldNames[idx_group].resize( 1 );	×
226
227	_outputFields[idx_group][0].resize( _nCells );	×
228	_outputFieldNames[idx_group][0] = "radiation flux density";	×
229	break;	×
230
231	case MEDICAL:	×
232	_outputFields[idx_group].resize( 2 );	×
233	_outputFieldNames[idx_group].resize( 2 );	×
234
235	// Dose
236	_outputFields[idx_group][0].resize( _nCells );	×
237	_outputFieldNames[idx_group][0] = "dose";	×
238	// Normalized Dose
239	_outputFields[idx_group][1].resize( _nCells );	×
240	_outputFieldNames[idx_group][1] = "normalized dose";	×
241	break;	×
242
243	default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;	×
244	}
245	}
246	}
247
248	void CSDSNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {	×
249	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();	×
250	double maxDose;
251	if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) \|\|	×
NEW 252	( idx_pseudoTime == _nIter - 1 ) /* need sol at last iteration */ ) {	×
253
254	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {	×
255	switch( _settings->GetVolumeOutput()[idx_group] ) {	×
256	case MINIMAL:	×
257	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
NEW 258	_outputFields[idx_group][0][idx_cell] = _scalarFluxNew[idx_cell];	×
259	}
260	break;	×
261
262	case MEDICAL:	×
263	// Compute Dose
264	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
265	_outputFields[idx_group][0][idx_cell] = _dose[idx_cell];	×
266	}
267	// Compute normalized dose
268	_outputFields[idx_group][1] = _outputFields[idx_group][0];	×
269
270	maxDose = *std::max_element( _outputFields[idx_group][0].begin(), _outputFields[idx_group][0].end() );	×
271
272	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
273	_outputFields[idx_group][1][idx_cell] /= maxDose;	×
274	}
275	break;	×
276
277	default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;	×
278	}
279	}
280	}
281	if( idx_pseudoTime == _nEnergies - 2 ) {	×
282	std::ofstream out( _settings->GetOutputFile().append( ".txt" ) );	×
283	unsigned nx = _settings->GetNCells();	×
284
285	for( unsigned j = 0; j < nx; ++j ) {	×
286	out << _cellMidPoints[j][0] << " " << _cellMidPoints[j][1] << " " << _dose[j] << std::endl;	×
287	}
288	out.close();	×
289	}
290	}

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