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

63.51

/src/problems/problembase.cpp

#include "problems/problembase.hpp"

#include "common/config.hpp"
#include "common/mesh.hpp"
#include "problems/aircavity1d.hpp"
#include "problems/checkerboard.hpp"
#include "problems/halflattice.hpp"
#include "problems/hohlraum.hpp"
#include "problems/lattice.hpp"
#include "problems/linesource.hpp"
#include "problems/meltingcube.hpp"
#include "problems/phantomimage.hpp"
#include "problems/quarterhohlraum.hpp"
#include "problems/radiationctimage.hpp"
#include "problems/starmapvalidation.hpp"
#include "problems/symmetrichohlraum.hpp"
#include "quadratures/quadraturebase.hpp"
#include "toolboxes/errormessages.hpp"

ProblemBase::ProblemBase( Config* settings, Mesh* mesh, QuadratureBase* quad ) {
    _settings = settings;
    _mesh     = mesh;
    _quad     = quad;

    // initialize QOI helper variables
    _curMaxOrdinateOutflow   = 0.0;
    _curScalarOutflow        = 0.0;
    _totalScalarOutflow      = 0.0;
    _mass                    = 0.0;
    _changeRateFlux          = 0.0;
    _dummyProbeMoments       = VectorVector( 4, Vector( 3, 0.0 ) );
    _dummyProbeValsGreenLine = std::vector<double>( _settings->GetNumProbingCellsLineHohlraum(), 0.0 );
    SetGhostCells();
}

ProblemBase::~ProblemBase() {}

ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh, QuadratureBase* quad ) {
    // Choose problem type
    switch( settings->GetProblemName() ) {
        case PROBLEM_Linesource: {
            if( settings->GetIsMomentSolver() )
                return new LineSource_Moment( settings, mesh, quad );
            else
                return new LineSource_SN( settings, mesh, quad );
        } break;
        case PROBLEM_Linesource1D: {
            if( settings->GetIsMomentSolver() )
                return new LineSource_Moment_1D( settings, mesh, quad );
            else
                return new LineSource_SN_1D( settings, mesh, quad );
        } break;
        case PROBLEM_Checkerboard: {
            if( settings->GetIsMomentSolver() )
                return new Checkerboard_Moment( settings, mesh, quad );
            else
                return new Checkerboard_SN( settings, mesh, quad );
        } break;
        case PROBLEM_Checkerboard1D: {
            if( settings->GetIsMomentSolver() )
                return new Checkerboard_Moment_1D( settings, mesh, quad );
            else
                return new Checkerboard_SN_1D( settings, mesh, quad );
        } break;
        case PROBLEM_Aircavity1D: {
            if( settings->GetIsMomentSolver() )
                return new AirCavity1D_Moment( settings, mesh, quad );
            else
                return new AirCavity1D( settings, mesh, quad );
        } break;
        case PROBLEM_StarmapValidation: {
            if( settings->GetIsMomentSolver() )
                return new StarMapValidation_Moment( settings, mesh, quad );
            else
                return new StarMapValidation_SN( settings, mesh, quad );
        } break;
        case PROBLEM_Phantomimage: return new PhantomImage( settings, mesh, quad );
        case PROBLEM_RadiationCT: {
            if( settings->GetIsMomentSolver() )
                return new RadiationCTImage_Moment( settings, mesh, quad );
            else
                return new RadiationCTImage( settings, mesh, quad );
        } break;
        case PROBLEM_Meltingcube: {
            if( settings->GetIsMomentSolver() )
                return new MeltingCube_Moment( settings, mesh, quad );
            else
                return new MeltingCube_SN( settings, mesh, quad );
        } break;
        case PROBLEM_Meltingcube1D: {
            if( settings->GetIsMomentSolver() )
                return new MeltingCube_Moment_1D( settings, mesh, quad );
            else
                return new MeltingCube_SN_1D( settings, mesh, quad );
        } break;
        case PROBLEM_Hohlraum: {
            if( settings->GetIsMomentSolver() )
                return new Hohlraum_Moment( settings, mesh, quad );
            else
                return new Hohlraum( settings, mesh, quad );
        } break;
        case PROBLEM_SymmetricHohlraum: {
            if( settings->GetIsMomentSolver() )
                return new SymmetricHohlraum_Moment( settings, mesh, quad );
            else
                return new SymmetricHohlraum( settings, mesh, quad );
        } break;
        case PROBLEM_QuarterHohlraum: {
            if( settings->GetIsMomentSolver() )
                return new QuarterHohlraum_Moment( settings, mesh, quad );
            else
                return new QuarterHohlraum( settings, mesh, quad );
        } break;
        case PROBLEM_Lattice: {
            if( settings->GetIsMomentSolver() )
                return new Lattice_Moment( settings, mesh, quad );
            else
                return new Lattice_SN( settings, mesh, quad );
        } break;
        case PROBLEM_HalfLattice: {
            if( settings->GetIsMomentSolver() )
                return new HalfLattice_Moment( settings, mesh, quad );
            else
                return new HalfLattice_SN( settings, mesh, quad );
        } break;

        default: ErrorMessages::Error( "No valid physical problem chosen. Please check your config file", CURRENT_FUNCTION ); return nullptr;
    }
}

// Default densities = 1
std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

// Legacy code: Scattering crossection loaded from database ENDF with physics
// class -> later overwritten with ICRU data
VectorVector ProblemBase::GetScatteringXSE( const Vector& /*energies*/, const Vector& /*angles*/ ) {
    ErrorMessages::Error( "Not yet implemented", CURRENT_FUNCTION );
    return VectorVector( 1, Vector( 1, 0 ) );
}

// Stopping powers from phyics class or default = -1
Vector ProblemBase::GetStoppingPower( const Vector& /* energies */ ) {
    ErrorMessages::Error( "Not yet implemented", CURRENT_FUNCTION );
    return Vector( 1, -1.0 );
}

void ProblemBase::SetGhostCells() {
    // Loop over all cells. If its a Dirichlet boundary, add cell to dict with
    // {cell_idx, boundary_value}
    auto cellBoundaries = _mesh->GetBoundaryTypes();
    std::map<int, Vector> ghostCellMap;

    Vector dummyGhostCell( 1, 0.0 );
    // #pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
        if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN || cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {
            // TODO: Refactor Boundary Conditions: We only have Ghost Cells with
            // Dirichlet conditions right now
            // #pragma omp critical
            { ghostCellMap.insert( { idx_cell, dummyGhostCell } ); }
        }
    }
    _ghostCells = ghostCellMap;
}

const Vector& ProblemBase::GetGhostCellValue( int /*idx_cell*/, const Vector& cell_sol ) { return cell_sol; }

void ProblemBase::ComputeMaxOrdinatewiseOutflow( const VectorVector& solution ) {
    if( _settings->GetSolverName() == SN_SOLVER || _settings->GetSolverName() == CSD_SN_SOLVER ) {
        double currOrdinatewiseOutflow = 0.0;
        double transportDirection      = 0.0;

        auto nCells   = _mesh->GetNumCells();
        auto cellMids = _mesh->GetCellMidPoints();
        auto areas    = _mesh->GetCellAreas();
        auto neigbors = _mesh->GetNeighbours();
        auto normals  = _mesh->GetNormals();

        auto quadPoints = _quad->GetPoints();
        auto weights    = _quad->GetWeights();
        auto nq         = _quad->GetNq();

        // Iterate over boundaries
        for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
            int idx_cell = it->first;    // Get Boundary cell index
            for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
                // Find face that points outward
                if( neigbors[idx_cell][idx_nbr] == nCells ) {
                    // Iterate over transport directions
                    for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {
                        transportDirection =
                            normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];
                        // Find outward facing transport directions
                        if( transportDirection > 0.0 ) {

                            currOrdinatewiseOutflow = transportDirection / norm( normals[idx_cell][idx_nbr] ) * solution[idx_cell][idx_quad];

                            if( currOrdinatewiseOutflow > _curMaxOrdinateOutflow ) _curMaxOrdinateOutflow = currOrdinatewiseOutflow;
                        }
                    }
                }
            }
        }
    }
    // TODO define alternative for Moment solvers
}

void ProblemBase::ComputeCurrentOutflow( const VectorVector& solution ) {
    if( _settings->GetSolverName() == SN_SOLVER || _settings->GetSolverName() == CSD_SN_SOLVER ) {

        _curScalarOutflow         = 0.0;
        double transportDirection = 0.0;

        auto nCells   = _mesh->GetNumCells();
        auto cellMids = _mesh->GetCellMidPoints();
        auto areas    = _mesh->GetCellAreas();
        auto neigbors = _mesh->GetNeighbours();
        auto normals  = _mesh->GetNormals();

        auto quadPoints = _quad->GetPoints();
        auto weights    = _quad->GetWeights();
        auto nq         = _quad->GetNq();

        // Iterate over boundaries
        for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
            int idx_cell = it->first;    // Get Boundary cell index

            // Iterate over face cell faces

            for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
                // Find face that points outward
                if( neigbors[idx_cell][idx_nbr] == nCells ) {
                    // Iterate over transport directions
                    for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {
                        transportDirection =
                            normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];
                        // Find outward facing transport directions
                        if( transportDirection > 0.0 ) {
                            _curScalarOutflow += transportDirection * solution[idx_cell][idx_quad] * weights[idx_quad];    // Integrate flux
                        }
                    }
                }
            }
        }
    }
    // TODO define alternative for Moment solvers
}

void ProblemBase::ComputeTotalOutflow( double dT ) { _totalScalarOutflow += _curScalarOutflow * dT; }

void ProblemBase::ComputeMass( const Vector& scalarFlux ) {
    _mass = 0.0;

    auto areas      = _mesh->GetCellAreas();
    unsigned nCells = _mesh->GetNumCells();
#pragma omp parallel for default( shared ) reduction( + : _mass )
    for( unsigned idx_cell = 0; idx_cell < nCells; ++idx_cell ) {
        _mass += scalarFlux[idx_cell] * areas[idx_cell];
    }
}

void ProblemBase::ComputeChangeRateFlux( const Vector& scalarFlux, const Vector& scalarFluxNew ) {
    _changeRateFlux = blaze::l2Norm( scalarFluxNew - scalarFlux );
}

1	#include "problems/problembase.hpp"
2
3	#include "common/config.hpp"
4	#include "common/mesh.hpp"
5	#include "problems/aircavity1d.hpp"
6	#include "problems/checkerboard.hpp"
7	#include "problems/halflattice.hpp"
8	#include "problems/hohlraum.hpp"
9	#include "problems/lattice.hpp"
10	#include "problems/linesource.hpp"
11	#include "problems/meltingcube.hpp"
12	#include "problems/phantomimage.hpp"
13	#include "problems/quarterhohlraum.hpp"
14	#include "problems/radiationctimage.hpp"
15	#include "problems/starmapvalidation.hpp"
16	#include "problems/symmetrichohlraum.hpp"
17	#include "quadratures/quadraturebase.hpp"
18	#include "toolboxes/errormessages.hpp"
19
20	ProblemBase::ProblemBase( Config* settings, Mesh* mesh, QuadratureBase* quad ) {	22✔
21	_settings = settings;	22✔
22	_mesh = mesh;	22✔
23	_quad = quad;	22✔
24
25	// initialize QOI helper variables
26	_curMaxOrdinateOutflow = 0.0;	22✔
27	_curScalarOutflow = 0.0;	22✔
28	_totalScalarOutflow = 0.0;	22✔
29	_mass = 0.0;	22✔
30	_changeRateFlux = 0.0;	22✔
31	_dummyProbeMoments = VectorVector( 4, Vector( 3, 0.0 ) );	22✔
32	_dummyProbeValsGreenLine = std::vector<double>( _settings->GetNumProbingCellsLineHohlraum(), 0.0 );	22✔
33	SetGhostCells();	22✔
34	}	22✔
35
36	ProblemBase::~ProblemBase() {}	16✔
37		×
38	ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh, QuadratureBase* quad ) {	16✔
39	// Choose problem type
40	switch( settings->GetProblemName() ) {	22✔
41	case PROBLEM_Linesource: {
42	if( settings->GetIsMomentSolver() )	22✔
43	return new LineSource_Moment( settings, mesh, quad );	12✔
44	else	12✔
45	return new LineSource_SN( settings, mesh, quad );	8✔
46	} break;
47	case PROBLEM_Linesource1D: {	4✔
48	if( settings->GetIsMomentSolver() )
NEW 49	return new LineSource_Moment_1D( settings, mesh, quad );	×
50	else	×
NEW 51	return new LineSource_SN_1D( settings, mesh, quad );	×
52	} break;
53	case PROBLEM_Checkerboard: {	×
54	if( settings->GetIsMomentSolver() )
55	return new Checkerboard_Moment( settings, mesh, quad );	6✔
56	else	6✔
57	return new Checkerboard_SN( settings, mesh, quad );	4✔
58	} break;
59	case PROBLEM_Checkerboard1D: {	2✔
60	if( settings->GetIsMomentSolver() )
NEW 61	return new Checkerboard_Moment_1D( settings, mesh, quad );	×
62	else	×
NEW 63	return new Checkerboard_SN_1D( settings, mesh, quad );	×
64	} break;
65	case PROBLEM_Aircavity1D: {	×
66	if( settings->GetIsMomentSolver() )
NEW 67	return new AirCavity1D_Moment( settings, mesh, quad );	×
68	else	×
NEW 69	return new AirCavity1D( settings, mesh, quad );	×
70	} break;
71	case PROBLEM_StarmapValidation: {	×
72	if( settings->GetIsMomentSolver() )
73	return new StarMapValidation_Moment( settings, mesh, quad );	4✔
74	else	4✔
75	return new StarMapValidation_SN( settings, mesh, quad );	4✔
76	} break;
NEW 77	case PROBLEM_Phantomimage: return new PhantomImage( settings, mesh, quad );	×
78	case PROBLEM_RadiationCT: {
79	if( settings->GetIsMomentSolver() )	×
NEW 80	return new RadiationCTImage_Moment( settings, mesh, quad );	×
81	else	×
NEW 82	return new RadiationCTImage( settings, mesh, quad );	×
83	} break;
84	case PROBLEM_Meltingcube: {	×
85	if( settings->GetIsMomentSolver() )
NEW 86	return new MeltingCube_Moment( settings, mesh, quad );	×
87	else	×
NEW 88	return new MeltingCube_SN( settings, mesh, quad );	×
89	} break;
90	case PROBLEM_Meltingcube1D: {	×
91	if( settings->GetIsMomentSolver() )
NEW 92	return new MeltingCube_Moment_1D( settings, mesh, quad );	×
93	else	×
NEW 94	return new MeltingCube_SN_1D( settings, mesh, quad );	×
95	} break;
96	case PROBLEM_Hohlraum: {	×
97	if( settings->GetIsMomentSolver() )
NEW 98	return new Hohlraum_Moment( settings, mesh, quad );	×
NEW 99	else	×
NEW 100	return new Hohlraum( settings, mesh, quad );	×
101	} break;
NEW 102	case PROBLEM_SymmetricHohlraum: {	×
103	if( settings->GetIsMomentSolver() )
NEW 104	return new SymmetricHohlraum_Moment( settings, mesh, quad );	×
NEW 105	else	×
NEW 106	return new SymmetricHohlraum( settings, mesh, quad );	×
107	} break;
NEW 108	case PROBLEM_QuarterHohlraum: {	×
109	if( settings->GetIsMomentSolver() )
NEW 110	return new QuarterHohlraum_Moment( settings, mesh, quad );	×
NEW 111	else	×
NEW 112	return new QuarterHohlraum( settings, mesh, quad );	×
113	} break;
NEW 114	case PROBLEM_Lattice: {	×
115	if( settings->GetIsMomentSolver() )
NEW 116	return new Lattice_Moment( settings, mesh, quad );	×
NEW 117	else	×
NEW 118	return new Lattice_SN( settings, mesh, quad );	×
119	} break;
NEW 120	case PROBLEM_HalfLattice: {	×
121	if( settings->GetIsMomentSolver() )
NEW 122	return new HalfLattice_Moment( settings, mesh, quad );	×
123	else	×
NEW 124	return new HalfLattice_SN( settings, mesh, quad );	×
125	} break;
126		×
127	default: ErrorMessages::Error( "No valid physical problem chosen. Please check your config file", CURRENT_FUNCTION ); return nullptr;
128	}
129	}	×
130
131	// Default densities = 1
132	std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
133
134	// Legacy code: Scattering crossection loaded from database ENDF with physics	36✔
135	// class -> later overwritten with ICRU data
136	VectorVector ProblemBase::GetScatteringXSE( const Vector& /energies/, const Vector& /angles/ ) {
137	ErrorMessages::Error( "Not yet implemented", CURRENT_FUNCTION );
138	return VectorVector( 1, Vector( 1, 0 ) );	×
139	}	×
140		×
141	// Stopping powers from phyics class or default = -1
142	Vector ProblemBase::GetStoppingPower( const Vector& /* energies */ ) {
143	ErrorMessages::Error( "Not yet implemented", CURRENT_FUNCTION );
144	return Vector( 1, -1.0 );	×
145	}	×
NEW 146		×
147	void ProblemBase::SetGhostCells() {
148	// Loop over all cells. If its a Dirichlet boundary, add cell to dict with
149	// {cell_idx, boundary_value}	22✔
150	auto cellBoundaries = _mesh->GetBoundaryTypes();
151	std::map<int, Vector> ghostCellMap;
152		44✔
153	Vector dummyGhostCell( 1, 0.0 );	44✔
154	// #pragma omp parallel for
155	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {	44✔
156	if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN \|\| cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {
157	// TODO: Refactor Boundary Conditions: We only have Ghost Cells with	11,110✔
158	// Dirichlet conditions right now	11,088✔
159	// #pragma omp critical
160	{ ghostCellMap.insert( { idx_cell, dummyGhostCell } ); }
161	}
162	}	1,184✔
163	_ghostCells = ghostCellMap;
164	}
165		22✔
166	const Vector& ProblemBase::GetGhostCellValue( int /idx_cell/, const Vector& cell_sol ) { return cell_sol; }	22✔
167
168	void ProblemBase::ComputeMaxOrdinatewiseOutflow( const VectorVector& solution ) {	1,008✔
169	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
170	double currOrdinatewiseOutflow = 0.0;	66✔
171	double transportDirection = 0.0;	66✔
172		22✔
173	auto nCells = _mesh->GetNumCells();	22✔
174	auto cellMids = _mesh->GetCellMidPoints();
175	auto areas = _mesh->GetCellAreas();	22✔
176	auto neigbors = _mesh->GetNeighbours();	44✔
177	auto normals = _mesh->GetNormals();	44✔
178		44✔
179	auto quadPoints = _quad->GetPoints();	44✔
180	auto weights = _quad->GetWeights();
181	auto nq = _quad->GetNq();	44✔
182		44✔
183	// Iterate over boundaries	22✔
184	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
185	int idx_cell = it->first; // Get Boundary cell index
186	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {	1,030✔
187	// Find face that points outward	1,008✔
188	if( neigbors[idx_cell][idx_nbr] == nCells ) {	4,032✔
189	// Iterate over transport directions
190	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	3,024✔
191	transportDirection =
192	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	87,696✔
193	// Find outward facing transport directions	86,688✔
194	if( transportDirection > 0.0 ) {	86,688✔
195
196	currOrdinatewiseOutflow = transportDirection / norm( normals[idx_cell][idx_nbr] ) * solution[idx_cell][idx_quad];	86,688✔
197
198	if( currOrdinatewiseOutflow > _curMaxOrdinateOutflow ) _curMaxOrdinateOutflow = currOrdinatewiseOutflow;	37,296✔
199	}
200	}	37,296✔
201	}
202	}
203	}
204	}
205	// TODO define alternative for Moment solvers
206	}
207
208	void ProblemBase::ComputeCurrentOutflow( const VectorVector& solution ) {	66✔
209	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
210		66✔
211	_curScalarOutflow = 0.0;	66✔
212	double transportDirection = 0.0;
213		22✔
214	auto nCells = _mesh->GetNumCells();	22✔
215	auto cellMids = _mesh->GetCellMidPoints();
216	auto areas = _mesh->GetCellAreas();	22✔
217	auto neigbors = _mesh->GetNeighbours();	44✔
218	auto normals = _mesh->GetNormals();	44✔
219		44✔
220	auto quadPoints = _quad->GetPoints();	44✔
221	auto weights = _quad->GetWeights();
222	auto nq = _quad->GetNq();	44✔
223		44✔
224	// Iterate over boundaries	22✔
225	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
226	int idx_cell = it->first; // Get Boundary cell index
227		1,030✔
228	// Iterate over face cell faces	1,008✔
229
230	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
231	// Find face that points outward
232	if( neigbors[idx_cell][idx_nbr] == nCells ) {	4,032✔
233	// Iterate over transport directions
234	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	3,024✔
235	transportDirection =
236	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	87,696✔
237	// Find outward facing transport directions	86,688✔
238	if( transportDirection > 0.0 ) {	86,688✔
239	_curScalarOutflow += transportDirection * solution[idx_cell][idx_quad] * weights[idx_quad]; // Integrate flux
240	}	86,688✔
241	}	37,296✔
242	}
243	}
244	}
245	}
246	// TODO define alternative for Moment solvers
247	}
248
249	void ProblemBase::ComputeTotalOutflow( double dT ) { _totalScalarOutflow += _curScalarOutflow * dT; }	66✔
250
251	void ProblemBase::ComputeMass( const Vector& scalarFlux ) {	66✔
252	_mass = 0.0;
253		66✔
254	auto areas = _mesh->GetCellAreas();	66✔
255	unsigned nCells = _mesh->GetNumCells();
256	#pragma omp parallel for default( shared ) reduction( + : _mass )	66✔
257	for( unsigned idx_cell = 0; idx_cell < nCells; ++idx_cell ) {	66✔
258	_mass += scalarFlux[idx_cell] * areas[idx_cell];	66✔
259	}
260	}
261
262	void ProblemBase::ComputeChangeRateFlux( const Vector& scalarFlux, const Vector& scalarFluxNew ) {	66✔
263	_changeRateFlux = blaze::l2Norm( scalarFluxNew - scalarFlux );
264	}	66✔

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