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

#include "problems/hohlraum.hpp"
#include "common/config.hpp"
#include "common/io.hpp"
#include "common/mesh.hpp"
#include "quadratures/qgausslegendretensorized.hpp"
#include "quadratures/quadraturebase.hpp"
#include "solvers/csdpn_starmap_constants.hpp"
#include "toolboxes/errormessages.hpp"
#include "toolboxes/interpolation.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"
#include "velocitybasis/sphericalbase.hpp"
#include "velocitybasis/sphericalharmonics.hpp"

Hohlraum::Hohlraum( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
    _sigmaS = Vector( _mesh->GetNumCells(), 0.1 );    // white area default
    _sigmaT = Vector( _mesh->GetNumCells(), 0.1 );    // white area default

#pragma omp parallel for
    for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
        double x = _mesh->GetCellMidPoints()[idx_cell][0];
        double y = _mesh->GetCellMidPoints()[idx_cell][1];
        // red area
        if( x < 0.05 && y > 0.25 && y < 1.05 ) {
            _sigmaS[idx_cell] = 95.0;
            _sigmaT[idx_cell] = 100.0;
        }
        // green area 1
        if( x > 0.45 && x < 0.85 && y > 0.25 && y < 0.3 ) {
            _sigmaS[idx_cell] = 90.0;
            _sigmaT[idx_cell] = 100.0;
        }
        // green area 2
        if( x > 0.45 && x < 0.85 && y > 1.0 && y < 1.05 ) {
            _sigmaS[idx_cell] = 90.0;
            _sigmaT[idx_cell] = 100.0;
        }
        // green area 3
        if( x > 0.45 && x < 0.5 && y > 0.25 && y < 1.05 ) {
            _sigmaS[idx_cell] = 90.0;
            _sigmaT[idx_cell] = 100.0;
        }
        // black area
        if( x > 0.5 && x < 0.85 && y > 0.3 && y < 1.0 ) {
            _sigmaS[idx_cell] = 50.0;
            _sigmaT[idx_cell] = 100.0;
        }
        // blue area
        if( x > 1.25 || y < 0.05 || y > 1.25 ) {
            _sigmaS[idx_cell] = 0.0;
            _sigmaT[idx_cell] = 100.0;
        }
    }
}

Hohlraum::~Hohlraum() {}

std::vector<VectorVector> Hohlraum::GetExternalSource( const Vector& /*energies*/ ) {
    VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
    auto cellMids = _mesh->GetCellMidPoints();
#pragma omp parallel for
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        // isotropic source at lef boundary region
        if( cellMids[j][0] < 0.05 ) {
            Q[j] = 0.0;    //_settings->GetSourceMagnitude();
        }
    }
    return std::vector<VectorVector>( 1u, Q );
}

VectorVector Hohlraum::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    VectorVector cellMids = _mesh->GetCellMidPoints();

    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        // boundary condition: Source on left side
        if( cellMids[j][0] < 0.0 && ( cellMids[j][1] > 0.0 && cellMids[j][1] < 1.3 ) ) {    // test case uses ghost cells
            psi[j] = _settings->GetSourceMagnitude();
            _mesh->SetBoundaryType( j, DIRICHLET );
        }
        else {
            psi[j] = 1e-4;
        }
    }
    return psi;
}

VectorVector Hohlraum::GetScatteringXS( const Vector& /*energies*/ ) { return VectorVector( 1u, _sigmaS ); }

VectorVector Hohlraum::GetTotalXS( const Vector& /*energies*/ ) { return VectorVector( 1u, _sigmaT ); }

Hohlraum_Moment::Hohlraum_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : Hohlraum( settings, mesh, quad ) {}

Hohlraum_Moment::~Hohlraum_Moment() {}

std::vector<VectorVector> Hohlraum_Moment::GetExternalSource( const Vector& /*energies*/ ) {
    // In case of PN, spherical basis is per default SPHERICAL_HARMONICS

    double integrationFactor = ( 4 * M_PI );
    if( _settings->GetDim() == 2 ) {
        integrationFactor = M_PI;
    }
    SphericalBase* tempBase  = SphericalBase::Create( _settings );
    unsigned ntotalEquations = tempBase->GetBasisSize();

    VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) );    // zero could lead to problems?
    VectorVector cellMids = _mesh->GetCellMidPoints();

    Vector uIC( ntotalEquations, 0 );

    if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
        QuadratureBase* quad          = QuadratureBase::Create( _settings );
        VectorVector quadPointsSphere = quad->GetPointsSphere();
        Vector w                      = quad->GetWeights();

        double my, phi;
        VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );

        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
            my                = quadPointsSphere[idx_quad][0];
            phi               = quadPointsSphere[idx_quad][1];
            moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
        }
        // Integrate <1*m> to get factors for monomial basis in isotropic scattering
        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
            uIC += w[idx_quad] * moments[idx_quad];
        }
        delete quad;
    }
    double kinetic_density = 0.0;    //_settings->GetSourceMagnitude();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( cellMids[j][0] < 0.05 ) {
            if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
                Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor;    // Remember scaling
            }
            if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                Q[j][0] = kinetic_density / integrationFactor;    // first bassis function is 1/ ( 4 * M_PI )
            }
        }
    }
    delete tempBase;    // Only temporally needed
    return std::vector<VectorVector>( 1u, Q );
}

VectorVector Hohlraum_Moment::SetupIC() {
    double integrationFactor = ( 4 * M_PI );
    if( _settings->GetDim() == 2 ) {
        integrationFactor = M_PI;
    }
    // In case of PN, spherical basis is per default SPHERICAL_HARMONICS
    SphericalBase* tempBase  = SphericalBase::Create( _settings );
    unsigned ntotalEquations = tempBase->GetBasisSize();

    VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );    // zero could lead to problems?
    VectorVector cellMids = _mesh->GetCellMidPoints();

    Vector tempIC( ntotalEquations, 0 );

    if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
        QuadratureBase* quad          = QuadratureBase::Create( _settings );
        VectorVector quadPointsSphere = quad->GetPointsSphere();
        Vector w                      = quad->GetWeights();

        double my, phi;
        VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );

        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
            my                = quadPointsSphere[idx_quad][0];
            phi               = quadPointsSphere[idx_quad][1];
            moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
        }
        // Integrate <1*m> to get factors for monomial basis in isotropic scattering
        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
            tempIC += w[idx_quad] * moments[idx_quad];
        }
        delete quad;
    }
    // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
    double kinetic_density = 1e-4;
    // std::vector<BOUNDARY_TYPE> _boundaryCells;
    for( unsigned j = 0; j < cellMids.size(); ++j ) {

        // boundary condition: Source on left side
        if( cellMids[j][0] < 0.0 && ( cellMids[j][1] > 0.0 && cellMids[j][1] < 1.3 ) ) {    // test case uses ghost cells
            kinetic_density = _settings->GetSourceMagnitude();
            _mesh->SetBoundaryType( j, DIRICHLET );
        }
        else {
            kinetic_density = 1e-4;
        }

        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
            initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor;    // Remember scaling
        }
        if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
            initialSolution[j][0] = kinetic_density / integrationFactor;    // first bassis function is 1/ ( 4 * M_PI )
        }
    }
    delete tempBase;    // Only temporally needed
    return initialSolution;
}

1	#include "problems/hohlraum.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	Hohlraum::Hohlraum( 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	#pragma omp parallel for	×
19	for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
20	double x = _mesh->GetCellMidPoints()[idx_cell][0];
21	double y = _mesh->GetCellMidPoints()[idx_cell][1];
22	// red area
23	if( x < 0.05 && y > 0.25 && y < 1.05 ) {
24	_sigmaS[idx_cell] = 95.0;
25	_sigmaT[idx_cell] = 100.0;
26	}
27	// green area 1
28	if( x > 0.45 && x < 0.85 && y > 0.25 && y < 0.3 ) {
29	_sigmaS[idx_cell] = 90.0;
30	_sigmaT[idx_cell] = 100.0;
31	}
32	// green area 2
33	if( x > 0.45 && x < 0.85 && y > 1.0 && y < 1.05 ) {
34	_sigmaS[idx_cell] = 90.0;
35	_sigmaT[idx_cell] = 100.0;
36	}
37	// green area 3
38	if( x > 0.45 && x < 0.5 && y > 0.25 && y < 1.05 ) {
39	_sigmaS[idx_cell] = 90.0;
40	_sigmaT[idx_cell] = 100.0;
41	}
42	// black area
43	if( x > 0.5 && x < 0.85 && y > 0.3 && y < 1.0 ) {
44	_sigmaS[idx_cell] = 50.0;
45	_sigmaT[idx_cell] = 100.0;
46	}
47	// blue area
48	if( x > 1.25 \|\| y < 0.05 \|\| y > 1.25 ) {
49	_sigmaS[idx_cell] = 0.0;
50	_sigmaT[idx_cell] = 100.0;
51	}
52	}
53	}
54
55	Hohlraum::~Hohlraum() {}	×
56		×
NEW 57	std::vector<VectorVector> Hohlraum::GetExternalSource( const Vector& /energies/ ) {	×
58	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
59	auto cellMids = _mesh->GetCellMidPoints();	×
60	#pragma omp parallel for	×
61	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
62	// isotropic source at lef boundary region	×
63	if( cellMids[j][0] < 0.05 ) {
64	Q[j] = 0.0; //_settings->GetSourceMagnitude();
65	}
66	}
67	return std::vector<VectorVector>( 1u, Q );
68	}
69		×
70	VectorVector Hohlraum::SetupIC() {
71	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
72	VectorVector cellMids = _mesh->GetCellMidPoints();	×
73		×
74	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
75	// boundary condition: Source on left side
76	if( cellMids[j][0] < 0.0 && ( cellMids[j][1] > 0.0 && cellMids[j][1] < 1.3 ) ) { // test case uses ghost cells	×
77	psi[j] = _settings->GetSourceMagnitude();
78	_mesh->SetBoundaryType( j, DIRICHLET );	×
79	}	×
80	else {	×
81	psi[j] = 1e-4;
82	}
83	}	×
84	return psi;
85	}
86		×
87	VectorVector Hohlraum::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }
88
NEW 89	VectorVector Hohlraum::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaT ); }	×
90
NEW 91	Hohlraum_Moment::Hohlraum_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : Hohlraum( settings, mesh, quad ) {}	×
92
93	Hohlraum_Moment::~Hohlraum_Moment() {}	×
94
NEW 95	std::vector<VectorVector> Hohlraum_Moment::GetExternalSource( const Vector& /energies/ ) {	×
96	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
97		×
98	double integrationFactor = ( 4 * M_PI );
99	if( _settings->GetDim() == 2 ) {	×
100	integrationFactor = M_PI;
101	}
102	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
103	unsigned ntotalEquations = tempBase->GetBasisSize();	×
104		×
105	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
106	VectorVector cellMids = _mesh->GetCellMidPoints();	×
107		×
108	Vector uIC( ntotalEquations, 0 );
109		×
NEW 110	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
111	QuadratureBase* quad = QuadratureBase::Create( _settings );
112	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
113	Vector w = quad->GetWeights();
114		×
115	double my, phi;	×
116	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
117		×
118	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
119	my = quadPointsSphere[idx_quad][0];
120	phi = quadPointsSphere[idx_quad][1];	×
121	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
122	}	×
123	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
124	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
125	uIC += w[idx_quad] * moments[idx_quad];	×
126	}
127	delete quad;
128	}	×
129	double kinetic_density = 0.0; //_settings->GetSourceMagnitude();	×
130	for( unsigned j = 0; j < cellMids.size(); ++j ) {
131	if( cellMids[j][0] < 0.05 ) {	×
132	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
133	Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor; // Remember scaling	×
134	}	×
135	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
136	Q[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
137	}	×
138	}
139	}	×
140	delete tempBase; // Only temporally needed	×
141	return std::vector<VectorVector>( 1u, Q );
142	}
143
144	VectorVector Hohlraum_Moment::SetupIC() {	×
145	double integrationFactor = ( 4 * M_PI );	×
146	if( _settings->GetDim() == 2 ) {
147	integrationFactor = M_PI;
148	}	×
149	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
150	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
151	unsigned ntotalEquations = tempBase->GetBasisSize();	×
152
153	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
154	VectorVector cellMids = _mesh->GetCellMidPoints();	×
155		×
156	Vector tempIC( ntotalEquations, 0 );
157		×
NEW 158	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
159	QuadratureBase* quad = QuadratureBase::Create( _settings );
160	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
161	Vector w = quad->GetWeights();
162		×
163	double my, phi;	×
164	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
165		×
166	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
167	my = quadPointsSphere[idx_quad][0];
168	phi = quadPointsSphere[idx_quad][1];	×
169	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
170	}	×
171	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
172	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
173	tempIC += w[idx_quad] * moments[idx_quad];	×
174	}
175	delete quad;
176	}	×
177	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
178	double kinetic_density = 1e-4;
179	// std::vector<BOUNDARY_TYPE> _boundaryCells;	×
180	for( unsigned j = 0; j < cellMids.size(); ++j ) {
181
182	// boundary condition: Source on left side	×
183	if( cellMids[j][0] < 0.0 && ( cellMids[j][1] > 0.0 && cellMids[j][1] < 1.3 ) ) { // test case uses ghost cells
184	kinetic_density = _settings->GetSourceMagnitude();	×
185	_mesh->SetBoundaryType( j, DIRICHLET );
186	}
187	else {	×
188	kinetic_density = 1e-4;	×
189	}	×
190
191	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
192	initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor; // Remember scaling	×
193	}
194	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
195	initialSolution[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	×
196	}	×
197	}
198	delete tempBase; // Only temporally needed	×
199	return initialSolution;	×
200	}

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