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

1	#include "problems/quarterhohlraum.hpp"
2	#include "common/config.hpp"
3	#include "common/io.hpp"
4	#include "common/mesh.hpp"
5	#include "quadratures/qgausslegendretensorized.hpp"
6	#include "quadratures/quadraturebase.hpp"
7	#include "solvers/csdpn_starmap_constants.hpp"
8	#include "toolboxes/errormessages.hpp"
9	#include "toolboxes/interpolation.hpp"
10	#include "toolboxes/textprocessingtoolbox.hpp"
11	#include "velocitybasis/sphericalbase.hpp"
12	#include "velocitybasis/sphericalharmonics.hpp"
13
NEW 14	QuarterHohlraum::QuarterHohlraum( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {	×
NEW 15	_sigmaS = Vector( _mesh->GetNumCells(), 0.1 ); // white area default	×
NEW 16	_sigmaT = Vector( _mesh->GetNumCells(), 0.1 ); // white area default	×
17
18	// Geometry of the red barrier
NEW 19	_redRightTop = _settings->GetPosRedRightTopHohlraum();	×
NEW 20	_posRedRightBorder = _settings->GetPosRedRightBorderHohlraum();	×
21
22	// Geometry of the green capsule
NEW 23	_thicknessGreen = 0.05;	×
NEW 24	_cornerUpperLeftGreen = { 0., 0.4 - _thicknessGreen / 2.0 };	×
NEW 25	_cornerLowerLeftGreen = { 0., +_thicknessGreen / 2.0 };	×
NEW 26	_cornerUpperRightGreen = { 0.2 - _thicknessGreen / 2.0, 0.4 - _thicknessGreen / 2.0 };	×
NEW 27	_cornerLowerRightGreen = { 0.2 - _thicknessGreen / 2.0, 0. + _thicknessGreen / 2.0 };	×
28
NEW 29	_curAbsorptionHohlraumCenter = 0.0;	×
NEW 30	_curAbsorptionHohlraumVertical = 0.0;	×
NEW 31	_curAbsorptionHohlraumHorizontal = 0.0;	×
NEW 32	_totalAbsorptionHohlraumCenter = 0.0;	×
NEW 33	_totalAbsorptionHohlraumVertical = 0.0;	×
NEW 34	_totalAbsorptionHohlraumHorizontal = 0.0;	×
NEW 35	_varAbsorptionHohlraumGreen = 0.0;	×
NEW 36	_probingCells = {	×
NEW 37	_mesh->GetCellOfKoordinate( 0.4, 0. ),	×
NEW 38	_mesh->GetCellOfKoordinate( 0., 0.5 ),	×
39	};
NEW 40	_probingMoments = VectorVector( 2, Vector( 3, 0.0 ) );	×
NEW 41	_nProbingCellsLineGreen = _settings->GetNumProbingCellsLineHohlraum();	×
NEW 42	SetProbingCellsLineGreen();	×
NEW 43	_absorptionValsIntegrated = std::vector<double>( _nProbingCellsLineGreen, 0.0 );	×
NEW 44	_varAbsorptionValsIntegrated = std::vector<double>( _nProbingCellsLineGreen, 0.0 );	×
45
NEW 46	#pragma omp parallel for	×
47	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
48	// Assumption: Domain size is 1.3x1.3
49	double x = _mesh->GetCellMidPoints()[idx_cell][0];
50	double y = _mesh->GetCellMidPoints()[idx_cell][1];
51
52	// red area right
53	if( x > _posRedRightBorder && y < _redRightTop ) {
54	_sigmaS[idx_cell] = 95.0;
55	_sigmaT[idx_cell] = 100.0;
56	}
57	// green and blue area
58	if( x > -0.2 && x < 0.2 && y > -0.4 && y < 0.4 ) {
59	_sigmaS[idx_cell] = 90.0;
60	_sigmaT[idx_cell] = 100.0;
61	}
62	// blue checkered area (overwrites part of green n blue area)
63	if( x > -0.15 && x < 0.15 && y > -0.35 && y < 0.35 ) {
64	_sigmaS[idx_cell] = 0.0;
65	_sigmaT[idx_cell] = 100.0;
66	}
67	// black area (upper and lower boundary)
68	if( y > 0.6 \|\| y < -0.6 ) {
69	_sigmaS[idx_cell] = 50.0;
70	_sigmaT[idx_cell] = 100.0;
71	}
72	}
NEW 73	SetGhostCells();	×
74	}
75
NEW 76	QuarterHohlraum::~QuarterHohlraum() {}	×
NEW 77		×
NEW 78	std::vector<VectorVector> QuarterHohlraum::GetExternalSource( const Vector& /* energies */ ) {	×
79	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
NEW 80	return std::vector<VectorVector>( 1u, Q );	×
NEW 81	}	×
NEW 82		×
83	VectorVector QuarterHohlraum::SetupIC() {
84	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
NEW 85	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 86		×
NEW 87	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
88	psi[j] = 0.0; // zero initial condition
NEW 89	}	×
NEW 90	return psi;	×
91	}
NEW 92		×
93	void QuarterHohlraum::SetGhostCells() {
94	// Loop over all cells. If its a Dirichlet boundary, add cell to dict with {cell_idx, boundary_value}
NEW 95	auto cellBoundaries = _mesh->GetBoundaryTypes();	×
96	std::map<int, Vector> ghostCellMap;
NEW 97	std::map<int, bool> ghostCellReflMap;	×
NEW 98		×
NEW 99	double tol = 1e-12; // For distance to boundary	×
100
NEW 101	QuadratureBase* quad = QuadratureBase::Create( _settings );	×
102	VectorVector vq = quad->GetPoints();
NEW 103	unsigned nq = quad->GetNq();	×
NEW 104		×
NEW 105	if( _settings->GetQuadName() != QUAD_GaussLegendreTensorized2D ) {	×
106	ErrorMessages::Error( "This simplified test case only works with symmetric quadrature orders. Use QUAD_GAUSS_LEGENDRE_TENSORIZED_2D",
NEW 107	CURRENT_FUNCTION );	×
NEW 108	}	×
109	{ // Create the symmetry maps for the quadratures
110
111	for( unsigned idx_q = 0; idx_q < nq; idx_q++ ) {
112	for( unsigned idx_q2 = 0; idx_q2 < nq; idx_q2++ ) {
NEW 113	if( abs( vq[idx_q][0] + vq[idx_q2][0] ) + abs( vq[idx_q][1] - vq[idx_q2][1] ) < tol ) {	×
NEW 114	_quadratureYReflection[idx_q] = idx_q2;	×
NEW 115	break;	×
NEW 116	}	×
NEW 117	}	×
118	for( unsigned idx_q2 = 0; idx_q2 < nq; idx_q2++ ) {
119	if( abs( vq[idx_q][0] - vq[idx_q2][0] ) + abs( vq[idx_q][1] + vq[idx_q2][1] ) < tol ) {
NEW 120	_quadratureXReflection[idx_q] = idx_q2;	×
NEW 121	break;	×
NEW 122	}	×
NEW 123	}	×
124	}
125	}
126	if( _quadratureXReflection.size() != nq ) {
127	ErrorMessages::Error( "Problem with X symmetry of quadrature of this mesh", CURRENT_FUNCTION );
NEW 128	}	×
NEW 129	if( _quadratureYReflection.size() != nq ) {	×
130	ErrorMessages::Error( "Problem with Y symmetry of quadrature of this mesh", CURRENT_FUNCTION );
NEW 131	}	×
NEW 132		×
133	Vector right_inflow( nq, 0.0 );
134	Vector vertical_flow( nq, 0.0 );
NEW 135		×
NEW 136	for( unsigned idx_q = 0; idx_q < nq; idx_q++ ) {	×
137	if( vq[idx_q][0] < 0.0 ) right_inflow[idx_q] = 1.0;
NEW 138	}	×
NEW 139		×
140	auto nodes = _mesh->GetNodes();
141	#pragma omp parallel for
NEW 142	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {	×
NEW 143	if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN \|\| cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {	×
144	#pragma omp critical
145	{
146	auto localCellNodes = _mesh->GetCells()[idx_cell];
147
148	_ghostCellsReflectingX[idx_cell] = false;
149	_ghostCellsReflectingY[idx_cell] = false;
150	for( unsigned idx_node = 0; idx_node < _mesh->GetNumNodesPerCell(); idx_node++ ) { // Check if corner node is in this cell
151	if( abs( nodes[localCellNodes[idx_node]][0] ) < tol ) { // close to 0 => left boundary
152	_ghostCellsReflectingY[idx_cell] = true;
153	ghostCellMap.insert( { idx_cell, vertical_flow } );
154	break;
155	}
156	else if( abs( nodes[localCellNodes[idx_node]][0] ) > 0.65 - tol ) { // right boundary
157	ghostCellMap.insert( { idx_cell, right_inflow } );
158	break;
159	}
160	else if( abs( nodes[localCellNodes[idx_node]][1] ) < tol ) { // lower boundary
161	_ghostCellsReflectingX[idx_cell] = true;
162	ghostCellMap.insert( { idx_cell, vertical_flow } );
163
164	break;
165	}
166	else if( abs( nodes[localCellNodes[idx_node]][1] ) > 0.65 - tol ) { // upper boundary
167	ghostCellMap.insert( { idx_cell, vertical_flow } );
168	break;
169	}
170	else if( idx_node == _mesh->GetNumNodesPerCell() - 1 ) {
171	ErrorMessages::Error( " Problem with ghost cell setup and boundary of this mesh at cell " + std::to_string( idx_cell ) +
172	" with node coordinates " + std::to_string( _mesh->GetNodes()[localCellNodes[idx_node]][0] ) + "," +
173	std::to_string( _mesh->GetNodes()[localCellNodes[idx_node]][1] ),
174	CURRENT_FUNCTION );
175	}
176	}
177	}
178	}
179	}
180	_ghostCells = ghostCellMap;
181
NEW 182	delete quad;	×
183	}
NEW 184		×
185	const Vector& QuarterHohlraum::GetGhostCellValue( int idx_cell, const Vector& cell_sol ) {
186	if( _ghostCellsReflectingX[idx_cell] ) {
NEW 187	for( unsigned idx_sys = 0; idx_sys < cell_sol.size(); idx_sys++ ) {	×
NEW 188	_ghostCells[idx_cell][idx_sys] = cell_sol[_quadratureXReflection[idx_sys]];	×
NEW 189	}	×
NEW 190	}	×
191	else if( _ghostCellsReflectingY[idx_cell] ) {
192	for( unsigned idx_sys = 0; idx_sys < cell_sol.size(); idx_sys++ ) {
NEW 193	_ghostCells[idx_cell][idx_sys] = cell_sol[_quadratureYReflection[idx_sys]];	×
NEW 194	}	×
NEW 195	}	×
196	return _ghostCells[idx_cell];
197	}
NEW 198		×
199	VectorVector QuarterHohlraum::GetScatteringXS( const Vector& /* energies */ ) { return VectorVector( 1u, _sigmaS ); }
200
NEW 201	VectorVector QuarterHohlraum::GetTotalXS( const Vector& /* energies */ ) { return VectorVector( 1u, _sigmaT ); }	×
202
NEW 203	void QuarterHohlraum::ComputeCurrentAbsorptionHohlraum( const Vector& scalarFlux ) {	×
204	_curAbsorptionHohlraumCenter = 0.0; // Green and blue areas of symmetric hohlraum
NEW 205	_curAbsorptionHohlraumVertical = 0.0; // Red areas of symmetric hohlraum	×
NEW 206	_curAbsorptionHohlraumHorizontal = 0.0; // Black areas of symmetric hohlraum	×
NEW 207		×
NEW 208	unsigned nCells = _mesh->GetNumCells();	×
209	auto cellMids = _mesh->GetCellMidPoints();
NEW 210	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 211		×
NEW 212	#pragma omp parallel for default( shared ) \	×
213	reduction( + : _curAbsorptionHohlraumCenter, _curAbsorptionHohlraumVertical, _curAbsorptionHohlraumHorizontal )
NEW 214	for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {	×
NEW 215	double x = _mesh->GetCellMidPoints()[idx_cell][0];	×
216	double y = _mesh->GetCellMidPoints()[idx_cell][1];
217
218	if( x > -0.2 && x < 0.2 && y > -0.35 && y < 0.35 ) {
219	_curAbsorptionHohlraumCenter += scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * areas[idx_cell];
220	}
221	if( ( x < -0.6 && y > -0.4 && y < 0.4 ) \|\| ( x > 0.6 && y > -0.4 && y < 0.4 ) ) {
222	_curAbsorptionHohlraumVertical += scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * areas[idx_cell];
223	}
224	if( y > 0.6 \|\| y < -0.6 ) {
225	_curAbsorptionHohlraumHorizontal += scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * areas[idx_cell];
226	}
227	}
228	}
229
230	void QuarterHohlraum::ComputeTotalAbsorptionHohlraum( double dT ) {
231	_totalAbsorptionHohlraumCenter += _curAbsorptionHohlraumCenter * dT;
NEW 232	_totalAbsorptionHohlraumVertical += _curAbsorptionHohlraumVertical * dT;	×
NEW 233	_totalAbsorptionHohlraumHorizontal += _curAbsorptionHohlraumHorizontal * dT;	×
NEW 234	}	×
NEW 235		×
236	void QuarterHohlraum::ComputeVarAbsorptionGreen( const Vector& scalarFlux ) {
237	double a_g = 0.0;
NEW 238	_varAbsorptionHohlraumGreen = 0.0;	×
NEW 239		×
NEW 240	unsigned nCells = _mesh->GetNumCells();	×
241	auto cellMids = _mesh->GetCellMidPoints();
NEW 242	std::vector<double> areas = _mesh->GetCellAreas();	×
NEW 243		×
NEW 244	#pragma omp parallel default( shared ) reduction( + : a_g )	×
245	for( unsigned idx_cell = 0; idx_cell < nCells; ++idx_cell ) {
NEW 246	double x = cellMids[idx_cell][0];	×
247	double y = cellMids[idx_cell][1];
248	bool green1, green2, green3, green4;
249
250	green1 = x > -0.2 && x < -0.15 && y > -0.35 && y < 0.35; // green area 1 (lower boundary)
251	green2 = x > 0.15 && x < 0.2 && y > -0.35 && y < 0.35; // green area 2 (upper boundary)
252	green3 = x > -0.2 && x < 0.2 && y > -0.4 && y < -0.35; // green area 3 (left boundary)
253	green4 = x > -0.2 && x < 0.2 && y > 0.35 && y < 0.4; // green area 4 (right boundary)
254
255	if( green1 \|\| green2 \|\| green3 \|\| green4 ) {
256	a_g += ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) * scalarFlux[idx_cell] * areas[idx_cell];
257	}
258	}
259
260	#pragma omp parallel default( shared ) reduction( + : _varAbsorptionHohlraumGreen )
261	for( unsigned idx_cell = 0; idx_cell < nCells; ++idx_cell ) {
NEW 262	double x = cellMids[idx_cell][0];	×
263	double y = cellMids[idx_cell][1];
264	bool green1, green2, green3, green4;
265
266	green1 = x > -0.2 && x < -0.15 && y > -0.35 && y < 0.35; // green area 1 (lower boundary)
267	green2 = x > 0.15 && x < 0.2 && y > -0.35 && y < 0.35; // green area 2 (upper boundary)
268	green3 = x > -0.2 && x < 0.2 && y > -0.4 && y < -0.35; // green area 3 (left boundary)
269	green4 = x > -0.2 && x < 0.2 && y > 0.35 && y < 0.4; // green area 4 (right boundary)
270
271	if( green1 \|\| green2 \|\| green3 \|\| green4 ) {
272	_varAbsorptionHohlraumGreen += ( a_g - scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) ) *
273	( a_g - scalarFlux[idx_cell] * ( _sigmaT[idx_cell] - _sigmaS[idx_cell] ) ) * areas[idx_cell];
274	}
275	}
276	}
277
278	void QuarterHohlraum::ComputeCurrentProbeMoment( const VectorVector& solution ) {
279	const VectorVector& quadPoints = _quad->GetPoints();
NEW 280	const Vector& weights = _quad->GetWeights();	×
NEW 281	unsigned nq = _quad->GetNq();	×
NEW 282		×
NEW 283	for( unsigned idx_cell = 0; idx_cell < 2; idx_cell++ ) { // Loop over probing cells	×
284	_probingMoments[idx_cell][0] = blaze::dot( solution[_probingCells[idx_cell]], weights );
NEW 285	_probingMoments[idx_cell][1] = 0.0;	×
NEW 286	_probingMoments[idx_cell][2] = 0.0;	×
NEW 287		×
NEW 288	for( unsigned idx_quad = 0; idx_quad < nq; idx_quad++ ) {	×
289	_probingMoments[idx_cell][1] += quadPoints[idx_quad][0] * solution[_probingCells[idx_cell]][idx_quad] * weights[idx_quad];
NEW 290	_probingMoments[idx_cell][2] += quadPoints[idx_quad][2] * solution[_probingCells[idx_cell]][idx_quad] * weights[idx_quad];	×
NEW 291	}	×
NEW 292	}	×
293	}
294
295	void QuarterHohlraum::SetProbingCellsLineGreen() {
296
NEW 297	double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] );	×
298	double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] );
NEW 299		×
NEW 300	// double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen );	×
301
302	unsigned nHorizontalProbingCells =
303	(unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) );
NEW 304	unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells;	×
NEW 305		×
NEW 306	_probingCellsLineGreen = std::vector<unsigned>( _nProbingCellsLineGreen );	×
307
NEW 308	// printf( "here" );	×
309
310	// Sample points on each side of the rectangle
311	std::vector<unsigned> side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells );
312	std::vector<unsigned> side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells );
NEW 313		×
NEW 314	// printf( "here" );	×
315	// Combine the points from each side
316	_probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() );
317	_probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() );
NEW 318	}	×
NEW 319		×
320	void QuarterHohlraum::ComputeQOIsGreenProbingLine( const Vector& scalarFlux ) {
321
NEW 322	double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] );	×
323	double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] );
NEW 324		×
NEW 325	double dl = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen );	×
326	double area = dl * _thicknessGreen;
NEW 327	double a_g = 0;	×
NEW 328	double l_max = _nProbingCellsLineGreen * dl;	×
NEW 329		×
NEW 330	for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells	×
331	_absorptionValsIntegrated[i] =
NEW 332	( _sigmaT[_probingCellsLineGreen[i]] - _sigmaS[_probingCellsLineGreen[i]] ) * scalarFlux[_probingCellsLineGreen[i]] * area;	×
NEW 333	a_g += _absorptionValsIntegrated[i] / (double)_nProbingCellsLineGreen;	×
NEW 334	}	×
NEW 335	for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells	×
336	_varAbsorptionValsIntegrated[i] = dl / l_max * ( a_g - _absorptionValsIntegrated[i] ) * ( a_g - _absorptionValsIntegrated[i] );
NEW 337	}	×
NEW 338	}	×
339
340	std::vector<unsigned> QuarterHohlraum::linspace2D( const std::vector<double>& start, const std::vector<double>& end, unsigned num_points ) {
341	std::vector<unsigned> result;
NEW 342	result.resize( num_points );	×
NEW 343	double stepX = ( end[0] - start[0] ) / ( num_points - 1 );	×
NEW 344	double stepY = ( end[1] - start[1] ) / ( num_points - 1 );	×
NEW 345		×
NEW 346	for( unsigned i = 0; i < num_points; ++i ) {	×
347	double x = start[0] + i * stepX;
NEW 348	double y = start[1] + i * stepY;	×
NEW 349		×
NEW 350	result[i] = _mesh->GetCellOfKoordinate( x, y );	×
351	}
NEW 352		×
353	return result;
354	}
NEW 355		×
356	void QuarterHohlraum::ComputeCurrentOutflow( const VectorVector& solution ) {
357	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
NEW 358		×
NEW 359	_curScalarOutflow = 0.0;	×
360	double transportDirection = 0.0;
NEW 361		×
NEW 362	auto nCells = _mesh->GetNumCells();	×
363	auto cellMids = _mesh->GetCellMidPoints();
NEW 364	auto areas = _mesh->GetCellAreas();	×
NEW 365	auto neigbors = _mesh->GetNeighbours();	×
NEW 366	auto normals = _mesh->GetNormals();	×
NEW 367		×
NEW 368	auto quadPoints = _quad->GetPoints();	×
369	auto weights = _quad->GetWeights();
NEW 370	auto nq = _quad->GetNq();	×
NEW 371		×
NEW 372	// Iterate over boundaries	×
373	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
374	int idx_cell = it->first; // Get Boundary cell index
NEW 375		×
NEW 376	// Iterate over face cell faces	×
377	if( !_ghostCellsReflectingX[idx_cell] && !_ghostCellsReflectingY[idx_cell] ) {
378	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {
NEW 379	// Find face that points outward	×
NEW 380	if( neigbors[idx_cell][idx_nbr] == nCells ) {	×
381	// Iterate over transport directions
NEW 382	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	×
383	transportDirection =
NEW 384	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	×
NEW 385	// Find outward facing transport directions	×
NEW 386	if( transportDirection > 0.0 ) {	×
387	_curScalarOutflow += transportDirection * solution[idx_cell][idx_quad] * weights[idx_quad]; // Integrate flux
NEW 388	}	×
NEW 389	}	×
390	}
391	}
392	}
393	}
394	}
395	// TODO define alternative for Moment solvers
396	}
397
398	void QuarterHohlraum::ComputeMaxOrdinatewiseOutflow( const VectorVector& solution ) {
399	if( _settings->GetSolverName() == SN_SOLVER \|\| _settings->GetSolverName() == CSD_SN_SOLVER ) {
NEW 400	double currOrdinatewiseOutflow = 0.0;	×
NEW 401	double transportDirection = 0.0;	×
NEW 402		×
NEW 403	auto nCells = _mesh->GetNumCells();	×
404	auto cellMids = _mesh->GetCellMidPoints();
NEW 405	auto areas = _mesh->GetCellAreas();	×
NEW 406	auto neigbors = _mesh->GetNeighbours();	×
NEW 407	auto normals = _mesh->GetNormals();	×
NEW 408		×
NEW 409	auto quadPoints = _quad->GetPoints();	×
410	auto weights = _quad->GetWeights();
NEW 411	auto nq = _quad->GetNq();	×
NEW 412		×
NEW 413	// Iterate over boundaries	×
414	for( std::map<int, Vector>::iterator it = _ghostCells.begin(); it != _ghostCells.end(); ++it ) {
415	int idx_cell = it->first; // Get Boundary cell index
NEW 416	if( !_ghostCellsReflectingX[idx_cell] && !_ghostCellsReflectingY[idx_cell] ) {	×
NEW 417	for( unsigned idx_nbr = 0; idx_nbr < neigbors[idx_cell].size(); ++idx_nbr ) {	×
NEW 418	// Find face that points outward	×
NEW 419	if( neigbors[idx_cell][idx_nbr] == nCells ) {	×
420	// Iterate over transport directions
NEW 421	for( unsigned idx_quad = 0; idx_quad < nq; ++idx_quad ) {	×
422	transportDirection =
NEW 423	normals[idx_cell][idx_nbr][0] * quadPoints[idx_quad][0] + normals[idx_cell][idx_nbr][1] * quadPoints[idx_quad][1];	×
NEW 424	// Find outward facing transport directions	×
NEW 425	if( transportDirection > 0.0 ) {	×
426	currOrdinatewiseOutflow = transportDirection / norm( normals[idx_cell][idx_nbr] ) * solution[idx_cell][idx_quad];
NEW 427	if( currOrdinatewiseOutflow > _curMaxOrdinateOutflow ) _curMaxOrdinateOutflow = currOrdinatewiseOutflow;	×
NEW 428	}	×
NEW 429	}	×
430	}
431	}
432	}
433	}
434	}
435	// TODO define alternative for Moment solvers
436	}
437	// -------------- Moment Symmetric Hohlraum ---------------
438
439	QuarterHohlraum_Moment::QuarterHohlraum_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : QuarterHohlraum( settings, mesh, quad ) {}
440
NEW 441	QuarterHohlraum_Moment::~QuarterHohlraum_Moment() {}	×
442
NEW 443	std::vector<VectorVector> QuarterHohlraum_Moment::GetExternalSource( const Vector& /* energies */ ) {	×
NEW 444	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
NEW 445		×
446	double integrationFactor = ( 4 * M_PI );
NEW 447	if( _settings->GetDim() == 2 ) {	×
448	integrationFactor = M_PI;
449	}
NEW 450	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 451	unsigned ntotalEquations = tempBase->GetBasisSize();	×
NEW 452		×
453	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
NEW 454	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 455		×
456	Vector uIC( ntotalEquations, 0 );
NEW 457		×
NEW 458	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
459	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 460	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
461	Vector w = quad->GetWeights();
NEW 462		×
NEW 463	double my, phi;	×
NEW 464	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 465		×
466	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
467	my = quadPointsSphere[idx_quad][0];
NEW 468	phi = quadPointsSphere[idx_quad][1];	×
469	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 470	}	×
NEW 471	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 472	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 473	uIC += w[idx_quad] * moments[idx_quad];	×
474	}
475	delete quad;
NEW 476	}	×
NEW 477	double kinetic_density = 0.0; //_settings->GetSourceMagnitude();	×
478	for( unsigned j = 0; j < cellMids.size(); ++j ) {
NEW 479	if( cellMids[j][0] < 0.05 ) {	×
480	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
NEW 481	Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor; // Remember scaling	×
NEW 482	}	×
NEW 483	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
NEW 484	Q[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 485	}	×
486	}
NEW 487	}	×
NEW 488	delete tempBase; // Only temporally needed	×
489	return std::vector<VectorVector>( 1u, Q );
490	}
491
NEW 492	VectorVector QuarterHohlraum_Moment::SetupIC() {	×
NEW 493	double integrationFactor = ( 4 * M_PI );	×
494	if( _settings->GetDim() == 2 ) {
495	integrationFactor = M_PI;
NEW 496	}	×
NEW 497	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
NEW 498	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
NEW 499	unsigned ntotalEquations = tempBase->GetBasisSize();	×
500
501	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
NEW 502	VectorVector cellMids = _mesh->GetCellMidPoints();	×
NEW 503		×
504	Vector tempIC( ntotalEquations, 0 );
NEW 505		×
NEW 506	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
507	QuadratureBase* quad = QuadratureBase::Create( _settings );
NEW 508	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
509	Vector w = quad->GetWeights();
NEW 510		×
NEW 511	double my, phi;	×
NEW 512	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
NEW 513		×
514	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
515	my = quadPointsSphere[idx_quad][0];
NEW 516	phi = quadPointsSphere[idx_quad][1];	×
517	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
NEW 518	}	×
NEW 519	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
NEW 520	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
NEW 521	tempIC += w[idx_quad] * moments[idx_quad];	×
522	}
523	delete quad;
NEW 524	}	×
NEW 525	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
526	double kinetic_density = 1e-4;
NEW 527	// std::vector<BOUNDARY_TYPE> _boundaryCells;	×
528	for( unsigned j = 0; j < cellMids.size(); ++j ) {
529
NEW 530	// boundary condition: Source on left side	×
531	if( cellMids[j][0] < 0.0 && ( cellMids[j][1] > 0.0 && cellMids[j][1] < 1.3 ) ) { // test case uses ghost cells
NEW 532	kinetic_density = _settings->GetSourceMagnitude();	×
533	_mesh->SetBoundaryType( j, DIRICHLET );
534	}
NEW 535	else {	×
NEW 536	kinetic_density = 1e-4;	×
NEW 537	}	×
538
539	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
NEW 540	initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor; // Remember scaling	×
541	}
542	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
NEW 543	initialSolution[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
NEW 544	}	×
545	}
NEW 546	delete tempBase; // Only temporally needed	×
NEW 547	return initialSolution;	×
548	}
549
NEW 550	void QuarterHohlraum_Moment::ComputeCurrentProbeMoment( const VectorVector& solution ) {	×
NEW 551	for( unsigned idx_cell = 0; idx_cell < 4; idx_cell++ ) { // Loop over probing cells	×
552	_probingMoments[idx_cell][0] = solution[_probingCells[idx_cell]][0];
553	if( _probingMoments[idx_cell].size() > 1 ) {
NEW 554	_probingMoments[idx_cell][1] = solution[_probingCells[idx_cell]][1];	×
NEW 555	_probingMoments[idx_cell][2] = solution[_probingCells[idx_cell]][2];	×
NEW 556	}	×
NEW 557	}	×
NEW 558	}	×

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