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

70.67

/src/solvers/snsolver.cpp

#include "solvers/snsolver.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/quadraturebase.hpp"
#include "toolboxes/errormessages.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"

#include "spdlog/spdlog.h"

SNSolver::SNSolver( Config* settings ) : SolverBase( settings ) {
    _quadPoints = _quadrature->GetPoints();
    _weights    = _quadrature->GetWeights();

    ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );
    _scatteringKernel   = k->GetScatteringKernel();
    delete k;

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

void SNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
    // Slope Limiter computation
    if( _reconsOrder > 1 ) {
        _mesh->ComputeSlopes( _nq, _solDx, _solDy, _sol );
        _mesh->ComputeLimiter( _nq, _solDx, _solDy, _sol, _limiter );
    }
}

void SNSolver::ComputeScalarFlux() {
    // double firstMomentScaleFactor = 4 * M_PI;
    // if( _settings->GetProblemName() == PROBLEM_Aircavity1D || _settings->GetProblemName() == PROBLEM_Linesource1D ||
    //    _settings->GetProblemName() == PROBLEM_Checkerboard1D || _settings->GetProblemName() == PROBLEM_Meltingcube1D ) {
    //    // firstMomentScaleFactor = 2.0;
    //}
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        _scalarFluxNew[idx_cell] = blaze::dot( _sol[idx_cell], _weights );    // / firstMomentScaleFactor;
    }
}

void SNSolver::FluxUpdate() {
    if( _settings->GetProblemName() == PROBLEM_Aircavity1D || _settings->GetProblemName() == PROBLEM_Linesource1D ||
        _settings->GetProblemName() == PROBLEM_Checkerboard1D || _settings->GetProblemName() == PROBLEM_Meltingcube1D ) {
        FluxUpdatePseudo1D();
    }
    else {
        FluxUpdatePseudo2D();
    }
}

void SNSolver::FluxUpdatePseudo1D() {
// Loop over all spatial cells
#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        double solL;
        double solR;
        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // Loop over all ordinates
        for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
            // Reset temporary variable
            _solNew[idx_cell][idx_quad] = 0.0;
            // Loop over all neighbor cells (edges) of cell j and compute numerical
            // fluxes
            for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
                // store flux contribution on psiNew_sigmaS to save memory
                if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells )
                    _solNew[idx_cell][idx_quad] +=
                        _g->Flux1D( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[idx_cell][idx_quad], _normals[idx_cell][idx_nbr] );
                else {
                    unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr];    // global idx of neighbor cell

                    switch( _reconsOrder ) {
                        // first order solver
                        case 1:
                            _solNew[idx_cell][idx_quad] +=
                                _g->Flux1D( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[nbr_glob][idx_quad], _normals[idx_cell][idx_nbr] );
                            break;
                        // second order solver
                        case 2:
                            // left status of interface
                            solL = _sol[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 = _sol[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] ) );
                            // flux evaluation
                            _solNew[idx_cell][idx_quad] += _g->Flux1D( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
                            break;
                            // higher order solver
                        default: ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION ); break;
                    }
                }
            }
        }
    }
}

void SNSolver::FluxUpdatePseudo2D() {
    if( _reconsOrder == 1 ) {
        // 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;
            // Reset temporary variable
            _solNew[idx_cell] *= 0.0;    // blaze op

            // Loop over all neighbor cells (edges) of cell j and compute numerical
            // fluxes
            for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
                // store flux contribution on psiNew_sigmaS to save memory
                if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
                    _g->Flux( _quadPoints,
                              _sol[idx_cell],
                              _problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
                              _solNew[idx_cell],
                              _normals[idx_cell][idx_nbr],
                              _nq );
                }
                else {
                    // first order solver
                    _g->Flux( _quadPoints, _sol[idx_cell], _sol[_neighbors[idx_cell][idx_nbr]], _solNew[idx_cell], _normals[idx_cell][idx_nbr], _nq );
                }
            }
        }
    }
    else if( _reconsOrder == 2 ) {
        // Loop over all spatial cells
#pragma omp parallel for
        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
            Vector solL;
            Vector solR;
            // Dirichlet cells stay at IC, farfield assumption
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
            // Loop over all ordinates
            _solNew[idx_cell] *= 0.0;    // blaze operation
                                         // Loop over all neighbor cells (edges) of cell j and compute numerical
                                         // fluxes
            for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
                // store flux contribution on psiNew_sigmaS to save memory
                if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
                    _g->Flux( _quadPoints,
                              _sol[idx_cell],
                              _problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
                              _solNew[idx_cell],
                              _normals[idx_cell][idx_nbr],
                              _nq );
                }
                else {
                    unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr];    // global idx of neighbor cell

                    // left status of interface
                    solL = _sol[idx_cell] +
                           _limiter[idx_cell] * ( _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
                                                  _solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );

                    solR = _sol[nbr_glob] +
                           _limiter[nbr_glob] * ( _solDx[nbr_glob] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
                                                  _solDy[nbr_glob] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );

                    // flux evaluation
                    _g->Flux( _quadPoints, solL, solR, _solNew[idx_cell], _normals[idx_cell][idx_nbr], _nq );
                }
            }
        }
    }
}

void SNSolver::FVMUpdate( unsigned idx_iter ) {
    if( _Q.size() == 1u && _sigmaT.size() == 1u && _sigmaS.size() == 1u ) {
#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;
            // loop over all ordinates
            // for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
            // time update angular flux with numerical flux and total scattering cross
            // section

            _solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] - _dT * _sigmaT[0][idx_cell] * _sol[idx_cell];
            //}
            // compute scattering effects
            _solNew[idx_cell] += _dT * _sigmaS[0][idx_cell] * _scatteringKernel * _sol[idx_cell];    // multiply scattering matrix with psi
            // Source Term
            // if( _Q[0][idx_cell].size() == 1u )    // isotropic source
            _solNew[idx_cell] += _dT * _Q[0][idx_cell][0];
            // else
            //     _solNew[idx_cell] += _dT * _Q[0][idx_cell];
        }
    }
    else {
#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;
            // loop over all ordinates
            // for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
            // time update angular flux with numerical flux and total scattering cross
            // section
            _solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] - _dT * _sigmaT[idx_iter][idx_cell] * _sol[idx_cell];
            // }
            // compute scattering effects
            _solNew[idx_cell] += _dT * _sigmaS[idx_iter][idx_cell] * _scatteringKernel * _sol[idx_cell];    // multiply scattering matrix with psi

            // Source Term
            if( _Q[0][idx_cell].size() == 1u )    // isotropic source
                _solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell][0];
            else
                _solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell];
        }
    }
}

void SNSolver::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 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 SN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

void SNSolver::WriteVolumeOutput( unsigned idx_iter ) {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
        ( idx_iter == _nEnergies - 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 ANALYTIC:
                    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 SN Solver!", CURRENT_FUNCTION ); break;
            }
        }
    }
}

1	#include "solvers/snsolver.hpp"
2
3	#include "common/config.hpp"
4	#include "common/io.hpp"
5	#include "common/mesh.hpp"
6	#include "fluxes/numericalflux.hpp"
7	#include "kernels/scatteringkernelbase.hpp"
8	#include "problems/problembase.hpp"
9	#include "quadratures/quadraturebase.hpp"
10	#include "toolboxes/errormessages.hpp"
11	#include "toolboxes/textprocessingtoolbox.hpp"
12
13	#include "spdlog/spdlog.h"
14
15	SNSolver::SNSolver( Config* settings ) : SolverBase( settings ) {	6✔
16	_quadPoints = _quadrature->GetPoints();	6✔
17	_weights = _quadrature->GetWeights();	6✔
18
19	ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );	6✔
20	_scatteringKernel = k->GetScatteringKernel();	6✔
21	delete k;	6✔
22
23	// Limiter variables
24	_solDx = VectorVector( _nCells, Vector( _nq, 0.0 ) );	6✔
25	_solDy = VectorVector( _nCells, Vector( _nq, 0.0 ) );	6✔
26	_limiter = VectorVector( _nCells, Vector( _nq, 0.0 ) );	6✔
27	}	6✔
28
29	void SNSolver::IterPreprocessing( unsigned /idx_iter/ ) {	22✔
30	// Slope Limiter computation
31	if( _reconsOrder > 1 ) {	22✔
32	_mesh->ComputeSlopes( _nq, _solDx, _solDy, _sol );	×
33	_mesh->ComputeLimiter( _nq, _solDx, _solDy, _sol, _limiter );	×
34	}
35	}	22✔
36
37	void SNSolver::ComputeScalarFlux() {	22✔
38	// double firstMomentScaleFactor = 4 * M_PI;
39	// if( _settings->GetProblemName() == PROBLEM_Aircavity1D \|\| _settings->GetProblemName() == PROBLEM_Linesource1D \|\|
40	// _settings->GetProblemName() == PROBLEM_Checkerboard1D \|\| _settings->GetProblemName() == PROBLEM_Meltingcube1D ) {
41	// // firstMomentScaleFactor = 2.0;
42	//}
43	#pragma omp parallel for	22✔
44	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
45	_scalarFluxNew[idx_cell] = blaze::dot( _sol[idx_cell], _weights ); // / firstMomentScaleFactor;
46	}
47	}	22✔
48
49	void SNSolver::FluxUpdate() {	22✔
50	if( _settings->GetProblemName() == PROBLEM_Aircavity1D \|\| _settings->GetProblemName() == PROBLEM_Linesource1D \|\|	66✔
51	_settings->GetProblemName() == PROBLEM_Checkerboard1D \|\| _settings->GetProblemName() == PROBLEM_Meltingcube1D ) {	66✔
52	FluxUpdatePseudo1D();	×
53	}
54	else {
55	FluxUpdatePseudo2D();	22✔
56	}
57	}	22✔
58
59	void SNSolver::FluxUpdatePseudo1D() {	×
60	// Loop over all spatial cells
61	#pragma omp parallel for	×
62	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
63	double solL;
64	double solR;
65	// Dirichlet cells stay at IC, farfield assumption
66	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
67	// Loop over all ordinates
68	for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
69	// Reset temporary variable
70	_solNew[idx_cell][idx_quad] = 0.0;
71	// Loop over all neighbor cells (edges) of cell j and compute numerical
72	// fluxes
73	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
74	// store flux contribution on psiNew_sigmaS to save memory
75	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells )
76	_solNew[idx_cell][idx_quad] +=
77	_g->Flux1D( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[idx_cell][idx_quad], _normals[idx_cell][idx_nbr] );
78	else {
79	unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
80
81	switch( _reconsOrder ) {
82	// first order solver
83	case 1:
84	_solNew[idx_cell][idx_quad] +=
85	_g->Flux1D( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[nbr_glob][idx_quad], _normals[idx_cell][idx_nbr] );
86	break;
87	// second order solver
88	case 2:
89	// left status of interface
90	solL = _sol[idx_cell][idx_quad] +
91	_limiter[idx_cell][idx_quad] *
92	( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
93	_solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
94	solR = _sol[nbr_glob][idx_quad] +
95	_limiter[nbr_glob][idx_quad] *
96	( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
97	_solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
98	// flux evaluation
99	_solNew[idx_cell][idx_quad] += _g->Flux1D( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
100	break;
101	// higher order solver
102	default: ErrorMessages::Error( "Reconstruction order not supported.", CURRENT_FUNCTION ); break;
103	}
104	}
105	}
106	}
107	}
108	}
109
110	void SNSolver::FluxUpdatePseudo2D() {	22✔
111	if( _reconsOrder == 1 ) {	22✔
112	// Loop over all spatial cells
113	#pragma omp parallel for	22✔
114	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
115	// Dirichlet cells stay at IC, farfield assumption
116	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
117	// Reset temporary variable
118	_solNew[idx_cell] *= 0.0; // blaze op
119
120	// Loop over all neighbor cells (edges) of cell j and compute numerical
121	// fluxes
122	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
123	// store flux contribution on psiNew_sigmaS to save memory
124	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
125	_g->Flux( _quadPoints,
126	_sol[idx_cell],
127	_problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
128	_solNew[idx_cell],
129	_normals[idx_cell][idx_nbr],
130	_nq );
131	}
132	else {
133	// first order solver
134	_g->Flux( _quadPoints, _sol[idx_cell], _sol[_neighbors[idx_cell][idx_nbr]], _solNew[idx_cell], _normals[idx_cell][idx_nbr], _nq );
135	}
136	}
137	}
138	}
NEW 139	else if( _reconsOrder == 2 ) {	×
140	// Loop over all spatial cells
NEW 141	#pragma omp parallel for	×
142	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
143	Vector solL;
144	Vector solR;
145	// Dirichlet cells stay at IC, farfield assumption
146	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
147	// Loop over all ordinates
148	_solNew[idx_cell] *= 0.0; // blaze operation
149	// Loop over all neighbor cells (edges) of cell j and compute numerical
150	// fluxes
151	for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
152	// store flux contribution on psiNew_sigmaS to save memory
153	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
154	_g->Flux( _quadPoints,
155	_sol[idx_cell],
156	_problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
157	_solNew[idx_cell],
158	_normals[idx_cell][idx_nbr],
159	_nq );
160	}
161	else {
162	unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
163
164	// left status of interface
165	solL = _sol[idx_cell] +
166	_limiter[idx_cell] * ( _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
167	_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
168
169	solR = _sol[nbr_glob] +
170	_limiter[nbr_glob] * ( _solDx[nbr_glob] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
171	_solDy[nbr_glob] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
172
173	// flux evaluation
174	_g->Flux( _quadPoints, solL, solR, _solNew[idx_cell], _normals[idx_cell][idx_nbr], _nq );
175	}
176	}
177	}
178	}
179	}	22✔
180
181	void SNSolver::FVMUpdate( unsigned idx_iter ) {	22✔
182	if( _Q.size() == 1u && _sigmaT.size() == 1u && _sigmaS.size() == 1u ) {	22✔
183	#pragma omp parallel for	22✔
184	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
185	// Dirichlet cells stay at IC, farfield assumption
186	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
187	// loop over all ordinates
188	// for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
189	// time update angular flux with numerical flux and total scattering cross
190	// section
191
192	_solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] - _dT * _sigmaT[0][idx_cell] * _sol[idx_cell];
193	//}
194	// compute scattering effects
195	_solNew[idx_cell] += _dT * _sigmaS[0][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi
196	// Source Term
197	// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
198	_solNew[idx_cell] += _dT * _Q[0][idx_cell][0];
199	// else
200	// _solNew[idx_cell] += _dT * _Q[0][idx_cell];
201	}
202	}
203	else {
NEW 204	#pragma omp parallel for	×
205	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
206	// Dirichlet cells stay at IC, farfield assumption
207	if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
208	// loop over all ordinates
209	// for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
210	// time update angular flux with numerical flux and total scattering cross
211	// section
212	_solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] - _dT * _sigmaT[idx_iter][idx_cell] * _sol[idx_cell];
213	// }
214	// compute scattering effects
215	_solNew[idx_cell] += _dT * _sigmaS[idx_iter][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi
216
217	// Source Term
218	if( _Q[0][idx_cell].size() == 1u ) // isotropic source
219	_solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell][0];
220	else
221	_solNew[idx_cell] += _dT * _Q[idx_iter][idx_cell];
222	}
223	}
224	}	22✔
225
226	void SNSolver::PrepareVolumeOutput() {	6✔
227	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();	6✔
228
229	_outputFieldNames.resize( nGroups );	6✔
230	_outputFields.resize( nGroups );	6✔
231
232	// Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
233	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {	12✔
234	// Prepare all Output Fields per group
235
236	// Different procedure, depending on the Group...
237	switch( _settings->GetVolumeOutput()[idx_group] ) {	6✔
238	case MINIMAL:	6✔
239	// Currently only one entry ==> rad flux
240	_outputFields[idx_group].resize( 1 );	6✔
241	_outputFieldNames[idx_group].resize( 1 );	6✔
242
243	_outputFields[idx_group][0].resize( _nCells );	6✔
244	_outputFieldNames[idx_group][0] = "radiation flux density";	6✔
245	break;	6✔
246
247	case ANALYTIC:	×
248	// one entry per cell
249	_outputFields[idx_group].resize( 1 );	×
250	_outputFieldNames[idx_group].resize( 1 );	×
251	_outputFields[idx_group][0].resize( _nCells );	×
252	_outputFieldNames[idx_group][0] = std::string( "analytic radiation flux density" );	×
253	break;	×
254
255	default: ErrorMessages::Error( "Volume Output Group not defined for SN Solver!", CURRENT_FUNCTION ); break;	×
256	}
257	}
258	}	6✔
259
260	void SNSolver::WriteVolumeOutput( unsigned idx_iter ) {	22✔
261	unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();	22✔
262	if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_iter % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) \|\|	44✔
263	( idx_iter == _nEnergies - 1 ) /* need sol at last iteration */ ) {	22✔
264	for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {	12✔
265	switch( _settings->GetVolumeOutput()[idx_group] ) {	6✔
266	case MINIMAL:	6✔
267	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	2,502✔
268	_outputFields[idx_group][0][idx_cell] = _scalarFluxNew[idx_cell];	2,496✔
269	}
270	break;	6✔
271	case ANALYTIC:	×
272	for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {	×
NEW 273	double time = idx_iter * _dT;	×
274	_outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(	×
NEW 275	_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_iter][idx_cell] );	×
276	}
277	break;	×
278
279	default: ErrorMessages::Error( "Volume Output Group not defined for SN Solver!", CURRENT_FUNCTION ); break;	×
280	}
281	}
282	}
283	}	22✔

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