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

1	#include "problems/halflattice.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	HalfLattice_SN::HalfLattice_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	×
NEW 23	auto log = spdlog::get( "event" );	×
NEW 24	log->info( "\| " );	×
NEW 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 {
NEW 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	HalfLattice_SN::~HalfLattice_SN() {}	×
NEW 59		×
NEW 60	VectorVector HalfLattice_SN::GetScatteringXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaS ); }	×
61
NEW 62	VectorVector HalfLattice_SN::GetTotalXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaT ); }	×
63
NEW 64	std::vector<VectorVector> HalfLattice_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 HalfLattice_SN::SetupIC() {
74	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
NEW 75	return psi;	×
NEW 76	}	×
NEW 77		×
78	bool HalfLattice_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	unsigned HalfLattice_SN::GetBlockID( const Vector& pos ) const {
95	double xy_corrector = 3.5;
NEW 96	int block_x = int( pos[0] + xy_corrector );	×
NEW 97	int block_y = int( pos[1] + xy_corrector );	×
NEW 98	return (unsigned)( block_y * 7 + block_x );	×
NEW 99	}	×
NEW 100		×
101	bool HalfLattice_SN::IsSource( const Vector& pos ) const {
102	// Check whether pos is part of source region
NEW 103	if( pos[0] >= 3 - 3.5 && pos[0] <= 4 - 3.5 && pos[1] >= 3 - 3.5 && pos[1] <= 4 - 3.5 )	×
104	return true;
NEW 105	else	×
NEW 106	return false;	×
107	}
NEW 108		×
109	void HalfLattice_SN::SetGhostCells() {
110	// Loop over all cells. If its a Dirichlet boundary, add cell to dict with {cell_idx, boundary_value}
NEW 111	auto cellBoundaries = _mesh->GetBoundaryTypes();	×
112	std::map<int, Vector> ghostCellMap;
NEW 113	std::map<int, bool> ghostCellReflMap;	×
NEW 114		×
NEW 115	double tol = 1e-12; // For distance to boundary	×
116
NEW 117	unsigned nGhostcells = 0;	×
118	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
NEW 119	if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN \|\| cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {	×
NEW 120	nGhostcells++;	×
NEW 121	}	×
NEW 122	}	×
123
124	QuadratureBase* quad = QuadratureBase::Create( _settings );
125	VectorVector vq = quad->GetPoints();
NEW 126	unsigned nq = quad->GetNq();	×
NEW 127		×
NEW 128	if( _settings->GetQuadName() != QUAD_GaussLegendreTensorized2D ) {	×
129	ErrorMessages::Error( "This simplified test case only works with symmetric quadrature orders. Use QUAD_GAUSS_LEGENDRE_TENSORIZED_2D",
NEW 130	CURRENT_FUNCTION );	×
NEW 131	}	×
132	{ // Create the symmetry maps for the quadratures
133
134	for( unsigned idx_q = 0; idx_q < nq; idx_q++ ) {
135	for( unsigned idx_q2 = 0; idx_q2 < nq; idx_q2++ ) {
NEW 136	if( abs( vq[idx_q][0] + vq[idx_q2][0] ) + abs( vq[idx_q][1] - vq[idx_q2][1] ) < tol ) {	×
NEW 137	_quadratureYReflection[idx_q] = idx_q2;	×
NEW 138	break;	×
NEW 139	}	×
NEW 140	}	×
141	}
142	}
143
144	if( _quadratureYReflection.size() != nq ) {
145	ErrorMessages::Error( "Problem with Y symmetry of quadrature of this mesh", CURRENT_FUNCTION );
NEW 146	}	×
NEW 147		×
148	Vector right_inflow( nq, 0.0 );
149	Vector vertical_flow( nq, 0.0 );
NEW 150		×
NEW 151	for( unsigned idx_q = 0; idx_q < nq; idx_q++ ) {	×
152	if( vq[idx_q][0] < 0.0 ) right_inflow[idx_q] = 1.0;
NEW 153	}	×
NEW 154		×
155	auto nodes = _mesh->GetNodes();
156
NEW 157	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {	×
158	if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN \|\| cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {
NEW 159		×
NEW 160	auto localCellNodes = _mesh->GetCells()[idx_cell];	×
161
NEW 162	_ghostCellsReflectingY[idx_cell] = false;	×
163	for( unsigned idx_node = 0; idx_node < _mesh->GetNumNodesPerCell(); idx_node++ ) { // Check if corner node is in this cell
NEW 164	if( abs( nodes[localCellNodes[idx_node]][1] ) < -3.5 + tol \|\| abs( nodes[localCellNodes[idx_node]][1] ) > 3.5 - tol ) {	×
NEW 165	// upper and lower boundary	×
NEW 166	ghostCellMap.insert( { idx_cell, vertical_flow } );	×
167	break;
NEW 168	}	×
NEW 169	else if( abs( nodes[localCellNodes[idx_node]][0] ) < tol ) { // close to 0 => left boundary	×
170	_ghostCellsReflectingY[idx_cell] = true;
NEW 171	ghostCellMap.insert( { idx_cell, vertical_flow } );	×
NEW 172	break;	×
NEW 173	}	×
NEW 174	else { // right boundary	×
175	ghostCellMap.insert( { idx_cell, vertical_flow } );
176	break;
NEW 177	}	×
NEW 178	}	×
179	}
180	}
181	_ghostCells = ghostCellMap;
182
NEW 183	delete quad;	×
184	}
NEW 185		×
186	const Vector& HalfLattice_SN::GetGhostCellValue( int idx_cell, const Vector& cell_sol ) {
187	if( _ghostCellsReflectingY[idx_cell] ) {
NEW 188	for( unsigned idx_sys = 0; idx_sys < cell_sol.size(); idx_sys++ ) {	×
NEW 189	_ghostCells[idx_cell][idx_sys] = cell_sol[_quadratureYReflection[idx_sys]];	×
NEW 190	}	×
NEW 191	}	×
192	return _ghostCells[idx_cell];
193	}
NEW 194		×
195	// QOI getter
196	double HalfLattice_SN::GetCurAbsorptionLattice() { return _curAbsorptionLattice; }
197	double HalfLattice_SN::GetTotalAbsorptionLattice() { return _totalAbsorptionLattice; }
NEW 198	double HalfLattice_SN::GetMaxAbsorptionLattice() { return _curMaxAbsorptionLattice; }	×
NEW 199		×
NEW 200	// QOI setter	×
201	void HalfLattice_SN::ComputeTotalAbsorptionLattice( double dT ) { _totalAbsorptionLattice += _curAbsorptionLattice * dT; }
202
NEW 203	void HalfLattice_SN::ComputeCurrentAbsorptionLattice( const Vector& scalarFlux ) {	×
204	_curAbsorptionLattice = 0.0;
NEW 205	unsigned nCells = _mesh->GetNumCells();	×
NEW 206	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 207	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 208		×
NEW 209	for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {	×
210	if( IsAbsorption( cellMids[idx_cell] ) ) {
NEW 211	_curAbsorptionLattice += scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * areas[idx_cell];	×
NEW 212	}	×
NEW 213	}	×
214	}
215	// TODO all absorption qois can be refactored in one function
216	void HalfLattice_SN::ComputeMaxAbsorptionLattice( const Vector& scalarFlux ) {
217	unsigned nCells = _mesh->GetNumCells();
NEW 218	auto cellMids = _mesh->GetCellMidPoints();	×
NEW 219	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 220		×
NEW 221	for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {	×
222	if( IsAbsorption( cellMids[idx_cell] ) ) {
NEW 223	if( _curMaxAbsorptionLattice < scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) )	×
NEW 224	_curMaxAbsorptionLattice = scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] );	×
NEW 225	}	×
NEW 226	}	×
227	}
228
229	void HalfLattice_SN::ComputeMaxOrdinatewiseOutflow( const VectorVector& solution ) {
230	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
NEW 231	double currOrdinatewiseOutflow = 0.0;	×
NEW 232	double transportDirection = 0.0;	×
NEW 233		×
NEW 234	auto nCells = _mesh->GetNumCells();	×
235	auto cellMids = _mesh->GetCellMidPoints();
NEW 236	auto areas = _mesh->GetCellAreas();	×
NEW 237	auto neigbors = _mesh->GetNeighbours();	×
NEW 238	auto normals = _mesh->GetNormals();	×
NEW 239		×
NEW 240	auto quadPoints = _quad->GetPoints();	×
241	auto weights = _quad->GetWeights();
NEW 242	auto nq = _quad->GetNq();	×
NEW 243		×
NEW 244	// Iterate over boundaries	×
245	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
246	int idx_cell = it->first; // Get Boundary cell index
NEW 247		×
NEW 248	if( !_ghostCellsReflectingY[idx_cell] ) { // Only work on non-reflecting boundaries	×
249	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
NEW 250	// Find face that points outward	×
NEW 251	if( neigbors[idx_cell][idx_nbr] == nCells ) {	×
252	// Iterate over transport directions
NEW 253	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	×
254	transportDirection =
NEW 255	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	×
NEW 256	// Find outward facing transport directions	×
NEW 257	if( transportDirection > 0.0 ) {	×
258
NEW 259	currOrdinatewiseOutflow = transportDirection / norm( normals[idx_cell][idx_nbr] ) * solution[idx_cell][idx_quad];	×
260
NEW 261	if( currOrdinatewiseOutflow > _curMaxOrdinateOutflow ) _curMaxOrdinateOutflow = currOrdinatewiseOutflow;	×
262	}
NEW 263	}	×
264	}
265	}
266	}
267	}
268	}
269	// TODO define alternative for Moment solvers
270	}
271
272	void HalfLattice_SN::ComputeCurrentOutflow( const VectorVector& solution ) {
273	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
NEW 274		×
NEW 275	_curScalarOutflow = 0.0;	×
276	double transportDirection = 0.0;
NEW 277		×
NEW 278	auto nCells = _mesh->GetNumCells();	×
279	auto cellMids = _mesh->GetCellMidPoints();
NEW 280	auto areas = _mesh->GetCellAreas();	×
NEW 281	auto neigbors = _mesh->GetNeighbours();	×
NEW 282	auto normals = _mesh->GetNormals();	×
NEW 283		×
NEW 284	auto quadPoints = _quad->GetPoints();	×
285	auto weights = _quad->GetWeights();
NEW 286	auto nq = _quad->GetNq();	×
NEW 287		×
NEW 288	// Iterate over boundaries	×
289	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
290	int idx_cell = it->first; // Get Boundary cell index
NEW 291		×
NEW 292	// Iterate over face cell faces	×
293	if( !_ghostCellsReflectingY[idx_cell] ) { // Only work on non-reflecting boundaries
294	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
NEW 295	// Find face that points outward	×
NEW 296	if( neigbors[idx_cell][idx_nbr] == nCells ) {	×
297	// Iterate over transport directions
NEW 298	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	×
299	transportDirection =
NEW 300	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	×
NEW 301	// Find outward facing transport directions	×
NEW 302	if( transportDirection > 0.0 ) {	×
303	_curScalarOutflow += transportDirection * solution[idx_cell][idx_quad] * weights[idx_quad]; // Integrate flux
NEW 304	}	×
NEW 305	}	×
306	}
307	}
308	}
309	}
310	}
311	// TODO define alternative for Moment solvers
312	}
313
314	// ---- Checkerboard Moments ----
315
316	// Constructor for checkerboard case with Pn
317	HalfLattice_Moment::HalfLattice_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
318
NEW 319	// Initialise crosssections = 1 (scattering)	×
320	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
321	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );
NEW 322		×
NEW 323	// for absorption regions change crosssections to all absorption	×
324	auto cellMids = _mesh->GetCellMidPoints();
325	for( unsigned j = 0; j < cellMids.size(); ++j ) {
NEW 326	if( isAbsorption( cellMids[j] ) ) {	×
NEW 327	_sigmaS[j] = 0.0;	×
NEW 328	_sigmaT[j] = 10.0;	×
NEW 329	}	×
NEW 330	}	×
331	}
332
333	HalfLattice_Moment::~HalfLattice_Moment() {}
334
NEW 335	VectorVector HalfLattice_Moment::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }	×
NEW 336		×
NEW 337	VectorVector HalfLattice_Moment::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaT ); }	×
338
NEW 339	std::vector<VectorVector> HalfLattice_Moment::GetExternalSource( const Vector& /energies/ ) {	×
340	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS
NEW 341		×
342	double integrationFactor = ( 4 * M_PI );
NEW 343	if( _settings->GetDim() == 2 ) {	×
344	integrationFactor = M_PI;
345	}
NEW 346	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 347	unsigned ntotalEquations = tempBase->GetBasisSize();	×
NEW 348		×
349	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
NEW 350	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 351		×
352	Vector uIC( ntotalEquations, 0 );
NEW 353		×
NEW 354	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
355	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 356	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
357	Vector w = quad->GetWeights();
NEW 358		×
NEW 359	double my, phi;	×
NEW 360	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 361		×
362	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
363	my = quadPointsSphere[idx_quad][0];
NEW 364	phi = quadPointsSphere[idx_quad][1];	×
365	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 366	}	×
NEW 367	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 368	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 369	uIC += w[idx_quad] * moments[idx_quad];	×
370	}
371	delete quad;
NEW 372	}	×
NEW 373	double kinetic_density = _settings->GetSourceMagnitude();	×
374	for( unsigned j = 0; j < cellMids.size(); ++j ) {
NEW 375	if( isSource( cellMids[j] ) ) {	×
376	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
NEW 377	Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor; // Remember scaling	×
NEW 378	}	×
NEW 379	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
NEW 380	Q[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 381	}	×
382	}
NEW 383	}	×
NEW 384	delete tempBase; // Only temporally needed	×
385	return std::vector<VectorVector>( 1u, Q );
386	}
387
NEW 388	VectorVector HalfLattice_Moment::SetupIC() {	×
NEW 389	double integrationFactor = ( 4 * M_PI );	×
390	if( _settings->GetDim() == 2 ) {
391	integrationFactor = M_PI;
NEW 392	}	×
NEW 393	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
NEW 394	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 395	unsigned ntotalEquations = tempBase->GetBasisSize();	×
396
397	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
NEW 398	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 399		×
400	Vector tempIC( ntotalEquations, 0 );
NEW 401		×
NEW 402	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
403	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 404	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
405	Vector w = quad->GetWeights();
NEW 406		×
NEW 407	double my, phi;	×
NEW 408	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 409		×
410	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
411	my = quadPointsSphere[idx_quad][0];
NEW 412	phi = quadPointsSphere[idx_quad][1];	×
413	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 414	}	×
NEW 415	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 416	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 417	tempIC += w[idx_quad] * moments[idx_quad];	×
418	}
419	delete quad;
NEW 420	}	×
NEW 421	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
422	double kinetic_density = 1e-4;
NEW 423	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
424	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
425	initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor; // Remember scaling
NEW 426	}	×
NEW 427	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
NEW 428	initialSolution[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 429	}	×
430	}
NEW 431	delete tempBase; // Only temporally needed	×
NEW 432	return initialSolution;	×
433	}
434
NEW 435	bool HalfLattice_Moment::isAbsorption( const Vector& pos ) const {	×
NEW 436	// Check whether pos is in absorption region	×
437	double xy_corrector = -3.5;
438	std::vector<double> lbounds{ 1 + xy_corrector, 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector };
NEW 439	std::vector<double> ubounds{ 2 + xy_corrector, 3 + xy_corrector, 4 + xy_corrector, 5 + xy_corrector, 6 + xy_corrector };	×
440	for( unsigned k = 0; k < lbounds.size(); ++k ) {
NEW 441	for( unsigned l = 0; l < lbounds.size(); ++l ) {	×
NEW 442	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	×
NEW 443	if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {	×
NEW 444	return true;	×
NEW 445	}	×
NEW 446	}	×
NEW 447	}	×
NEW 448	return false;	×
449	}
450
451	bool HalfLattice_Moment::isSource( const Vector& pos ) const {
NEW 452	// Check whether pos is in source region	×
453	if( pos[0] >= 3 - 3.5 && pos[0] <= 4 - 3.5 && pos[1] >= 3 - 3.5 && pos[1] <= 4 - 3.5 )
454	return true;
NEW 455	else	×
456	return false;
NEW 457	}	×

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