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/problems/lattice.cpp

1	#include "problems/lattice.hpp"
2	#include "common/config.hpp"
3	#include "common/mesh.hpp"
4	#include "quadratures/quadraturebase.hpp"
5	#include "toolboxes/textprocessingtoolbox.hpp"
6	#include "velocitybasis/sphericalbase.hpp"
7	#include "velocitybasis/sphericalharmonics.hpp"
8
9	// ---- Checkerboard Sn ----
10	// Constructor for Ckeckerboard case with Sn
NEW 11	Lattice_SN::Lattice_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {	×
12
13	// Initialise scattering crosssections to 1 and absorption cross sections to 0
NEW 14	_sigmaS = Vector( _mesh->GetNumCells(), 1. );	×
NEW 15	_sigmaT = Vector( _mesh->GetNumCells(), 1. );	×
16
17	// Initialize Quantities of interest
NEW 18	_curAbsorptionLattice = 0.0;	×
NEW 19	_curMaxAbsorptionLattice = 0.0;	×
NEW 20	_totalAbsorptionLattice = 0.0;	×
21
NEW 22	if( _settings->GetNLatticeAbsIndividual() == 49 && _settings->GetNLatticeScatterIndividual() == 49 ) { // Individual values set	×
23	// auto log = spdlog::get( "event" );
24	// log->info( "\| " );
25	// log->info( "\| Lattice test case WITH individual scattering and absorption values for each block. " );
26
NEW 27	auto cellMids = _mesh->GetCellMidPoints();	×
28
NEW 29	std::vector<double> scatteringValues = _settings->GetLatticeScatterIndividual();	×
NEW 30	std::vector<double> absorptionValues = _settings->GetLatticeAbsorptionIndividual();	×
31
NEW 32	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
NEW 33	_sigmaS[j] = scatteringValues[GetBlockID( cellMids[j] )];	×
NEW 34	_sigmaT[j] = absorptionValues[GetBlockID( cellMids[j] )] + scatteringValues[GetBlockID( cellMids[j] )];	×
35	}
36	}
37	else {
38	// auto log = spdlog::get( "event" );
39	// log->info( "\| " );
40	// log->info( "\| Lattice test case WITHOUT individual scattering and absorption values for each block. " );
41	// For absorption cells: set scattering XS to 0 and absorption to 10
NEW 42	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 43	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
NEW 44	if( IsAbsorption( cellMids[j] ) ) {	×
NEW 45	_sigmaS[j] = 0.0;	×
NEW 46	_sigmaT[j] = _settings->GetLatticeAbsBlue();	×
47	}
NEW 48	else if( !IsSource( cellMids[j] ) ) { // White block	×
NEW 49	_sigmaS[j] = _settings->GetLatticeScatterWhite();	×
NEW 50	_sigmaT[j] = _settings->GetLatticeScatterWhite();	×
51	}
52	}
53	}
54
NEW 55	SetGhostCells();	×
56	}
57
NEW 58	Lattice_SN::~Lattice_SN() {}	×
NEW 59		×
NEW 60	VectorVector Lattice_SN::GetScatteringXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaS ); }	×
61
NEW 62	VectorVector Lattice_SN::GetTotalXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaT ); }	×
63
NEW 64	std::vector<VectorVector> Lattice_SN::GetExternalSource( const Vector& /energies/ ) {	×
65	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
NEW 66	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 67	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
NEW 68	if( IsSource( cellMids[j] ) ) Q[j] = _settings->GetSourceMagnitude(); // isotropic source	×
NEW 69	}	×
NEW 70	return std::vector<VectorVector>( 1u, Q );	×
71	}
NEW 72		×
73	VectorVector Lattice_SN::SetupIC() {
74	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
NEW 75	return psi;	×
NEW 76	}	×
NEW 77		×
78	bool Lattice_SN::IsAbsorption( const Vector& pos ) const {
79	// Check whether pos is inside absorbing squares
NEW 80	double xy_corrector = -3.5;	×
81	std::vector<double> lbounds{ 1 + xy_corrector, 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector };
NEW 82	std::vector<double> ubounds{ 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector, 6 + xy_corrector };	×
NEW 83	for( unsigned k = 0; k < lbounds.size(); ++k ) {	×
NEW 84	for( unsigned l = 0; l < lbounds.size(); ++l ) {	×
NEW 85	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	×
NEW 86	if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {	×
NEW 87	return true;	×
NEW 88	}	×
NEW 89	}	×
90	}
91	return false;
92	}
NEW 93		×
94	bool Lattice_SN::IsAbsorption( double x, double y ) const {
95	// Check whether pos is inside absorbing squares
NEW 96	double xy_corrector = -3.5;	×
97	std::vector<double> lbounds{ 1 + xy_corrector, 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector };
NEW 98	std::vector<double> ubounds{ 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector, 6 + xy_corrector };	×
NEW 99	for( unsigned k = 0; k < lbounds.size(); ++k ) {	×
NEW 100	for( unsigned l = 0; l < lbounds.size(); ++l ) {	×
NEW 101	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	×
NEW 102	if( x >= lbounds[k] && x <= ubounds[k] && y >= lbounds[l] && y <= ubounds[l] ) {	×
NEW 103	return true;	×
NEW 104	}	×
NEW 105	}	×
106	}
107	return false;
108	}
NEW 109		×
110	unsigned Lattice_SN::GetBlockID( const Vector& pos ) const {
111	double xy_corrector = 3.5;
NEW 112	int block_x = int( pos[0] + xy_corrector );	×
NEW 113	int block_y = int( pos[1] + xy_corrector );	×
NEW 114	return (unsigned)( block_y * 7 + block_x );	×
NEW 115	}	×
NEW 116		×
117	bool Lattice_SN::IsSource( const Vector& pos ) const {
118	// Check whether pos is part of source region
NEW 119	if( pos[0] >= 3 - 3.5 && pos[0] <= 4 - 3.5 && pos[1] >= 3 - 3.5 && pos[1] <= 4 - 3.5 )	×
120	return true;
NEW 121	else	×
NEW 122	return false;	×
123	}
NEW 124		×
125	void Lattice_SN::SetGhostCells() {
126	// Loop over all cells. If its a Dirichlet boundary, add cell to dict with {cell_idx, boundary_value}
NEW 127	auto cellBoundaries = _mesh->GetBoundaryTypes();	×
128	std::map<int, Vector> ghostCellMap;
NEW 129		×
NEW 130	QuadratureBase* quad = QuadratureBase::Create( _settings );	×
131	VectorVector vq = quad->GetPoints();
NEW 132	unsigned nq = quad->GetNq();	×
NEW 133		×
NEW 134	Vector void_ghostcell( nq, 0.0 );	×
135
NEW 136	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {	×
137	if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN \|\| cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {
NEW 138	ghostCellMap[idx_cell] = void_ghostcell;	×
NEW 139	}	×
NEW 140	}	×
141	_ghostCells = ghostCellMap;
142
NEW 143	delete quad;	×
144	}
NEW 145		×
146	const Vector& Lattice_SN::GetGhostCellValue( int idx_cell, const Vector& /* cell_sol */ ) { return _ghostCells[idx_cell]; }
147
NEW 148	// QOI getter	×
149	double Lattice_SN::GetCurAbsorptionLattice() { return _curAbsorptionLattice; }
150	double Lattice_SN::GetTotalAbsorptionLattice() { return _totalAbsorptionLattice; }
NEW 151	double Lattice_SN::GetMaxAbsorptionLattice() { return _curMaxAbsorptionLattice; }	×
NEW 152	// QOI setter	×
NEW 153	void Lattice_SN::ComputeTotalAbsorptionLattice( double dT ) { _totalAbsorptionLattice += _curAbsorptionLattice * dT; }	×
154
NEW 155	void Lattice_SN::ComputeCurrentAbsorptionLattice( const Vector& scalarFlux ) {	×
156	_curAbsorptionLattice = 0.0;
NEW 157	unsigned nCells = _mesh->GetNumCells();	×
NEW 158	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 159	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 160		×
NEW 161	for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {	×
162	if( IsAbsorption( cellMids[idx_cell] ) ) {
NEW 163	_curAbsorptionLattice += scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * areas[idx_cell];	×
NEW 164	}	×
NEW 165	}	×
166	}
167	// TODO all absorption qois can be refactored in one function
168	void Lattice_SN::ComputeMaxAbsorptionLattice( const Vector& scalarFlux ) {
169	unsigned nCells = _mesh->GetNumCells();
NEW 170	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 171	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 172		×
NEW 173	for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {	×
174	if( IsAbsorption( cellMids[idx_cell] ) ) {
NEW 175	if( _curMaxAbsorptionLattice < scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) )	×
NEW 176	_curMaxAbsorptionLattice = scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] );	×
NEW 177	}	×
NEW 178	}	×
179	}
180
181	// ---- Checkerboard Moments ----
182
183	// Constructor for checkerboard case with Pn
184	Lattice_Moment::Lattice_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
185
NEW 186	// Initialise crosssections = 1 (scattering)	×
187	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
188	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );
NEW 189		×
NEW 190	// for absorption regions change crosssections to all absorption	×
191	auto cellMids = _mesh->GetCellMidPoints();
192	for( unsigned j = 0; j < cellMids.size(); ++j ) {
NEW 193	if( isAbsorption( cellMids[j] ) ) {	×
NEW 194	_sigmaS[j] = 0.0;	×
NEW 195	_sigmaT[j] = 10.0;	×
NEW 196	}	×
NEW 197	}	×
198	}
199
200	Lattice_Moment::~Lattice_Moment() {}
201
NEW 202	VectorVector Lattice_Moment::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }	×
NEW 203		×
NEW 204	VectorVector Lattice_Moment::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaT ); }	×
205
NEW 206	std::vector<VectorVector> Lattice_Moment::GetExternalSource( const Vector& /energies/ ) {	×
207	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS
NEW 208		×
209	double integrationFactor = ( 4 * M_PI );
NEW 210	if( _settings->GetDim() == 2 ) {	×
211	integrationFactor = M_PI;
212	}
NEW 213	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 214	unsigned ntotalEquations = tempBase->GetBasisSize();	×
NEW 215		×
216	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
NEW 217	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 218		×
219	Vector uIC( ntotalEquations, 0 );
NEW 220		×
NEW 221	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
222	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 223	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
224	Vector w = quad->GetWeights();
NEW 225		×
NEW 226	double my, phi;	×
NEW 227	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 228		×
229	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
230	my = quadPointsSphere[idx_quad][0];
NEW 231	phi = quadPointsSphere[idx_quad][1];	×
232	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 233	}	×
NEW 234	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 235	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 236	uIC += w[idx_quad] * moments[idx_quad];	×
237	}
238	delete quad;
NEW 239	}	×
NEW 240	double kinetic_density = _settings->GetSourceMagnitude();	×
241	for( unsigned j = 0; j < cellMids.size(); ++j ) {
NEW 242	if( isSource( cellMids[j] ) ) {	×
243	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
NEW 244	Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor; // Remember scaling	×
NEW 245	}	×
NEW 246	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
NEW 247	Q[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 248	}	×
249	}
NEW 250	}	×
NEW 251	delete tempBase; // Only temporally needed	×
252	return std::vector<VectorVector>( 1u, Q );
253	}
254
NEW 255	VectorVector Lattice_Moment::SetupIC() {	×
NEW 256	double integrationFactor = ( 4 * M_PI );	×
257	if( _settings->GetDim() == 2 ) {
258	integrationFactor = M_PI;
NEW 259	}	×
NEW 260	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
NEW 261	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 262	unsigned ntotalEquations = tempBase->GetBasisSize();	×
263
264	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
NEW 265	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 266		×
267	Vector tempIC( ntotalEquations, 0 );
NEW 268		×
NEW 269	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
270	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 271	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
272	Vector w = quad->GetWeights();
NEW 273		×
NEW 274	double my, phi;	×
NEW 275	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 276		×
277	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
278	my = quadPointsSphere[idx_quad][0];
NEW 279	phi = quadPointsSphere[idx_quad][1];	×
280	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 281	}	×
NEW 282	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 283	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 284	tempIC += w[idx_quad] * moments[idx_quad];	×
285	}
286	delete quad;
NEW 287	}	×
NEW 288	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
289	double kinetic_density = 1e-4;
NEW 290	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
291	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
292	initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor; // Remember scaling
NEW 293	}	×
NEW 294	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
NEW 295	initialSolution[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 296	}	×
297	}
NEW 298	delete tempBase; // Only temporally needed	×
NEW 299	return initialSolution;	×
300	}
301
NEW 302	bool Lattice_Moment::isAbsorption( const Vector& pos ) const {	×
NEW 303	// Check whether pos is in absorption region	×
304	double xy_corrector = -3.5;
305	std::vector<double> lbounds{ 1 + xy_corrector, 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector };
NEW 306	std::vector<double> ubounds{ 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector, 6 + xy_corrector };	×
307	for( unsigned k = 0; k < lbounds.size(); ++k ) {
NEW 308	for( unsigned l = 0; l < lbounds.size(); ++l ) {	×
NEW 309	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	×
NEW 310	if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {	×
NEW 311	return true;	×
NEW 312	}	×
NEW 313	}	×
NEW 314	}	×
NEW 315	return false;	×
316	}
317
318	bool Lattice_Moment::isSource( const Vector& pos ) const {
NEW 319	// Check whether pos is in source region	×
320	if( pos[0] >= 3 - 3.5 && pos[0] <= 4 - 3.5 && pos[1] >= 3 - 3.5 && pos[1] <= 4 - 3.5 )
321	return true;
NEW 322	else	×
323	return false;
NEW 324	}	×

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