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

#include "problems/meltingcube.hpp"
#include "common/config.hpp"
#include "common/mesh.hpp"
#include "quadratures/quadraturebase.hpp"
#include "velocitybasis/sphericalbase.hpp"
#include "velocitybasis/sphericalharmonics.hpp"
#include <complex>

// ---- Linesource ----

MeltingCube::MeltingCube( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
    _sigmaS = settings->GetSigmaS();
}

MeltingCube::~MeltingCube() {}

VectorVector MeltingCube::GetScatteringXS( const Vector& energies ) {
    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) );
}

VectorVector MeltingCube::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) ); }

std::vector<VectorVector> MeltingCube::GetExternalSource( const Vector& /*energies*/ ) {
    return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}

// ---- MeltingCube_SN ----

MeltingCube_SN::MeltingCube_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}

MeltingCube_SN::~MeltingCube_SN() {}

VectorVector MeltingCube_SN::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    auto cellMids          = _mesh->GetCellMidPoints();
    double kinetic_density = 0.0;
    double epsilon         = 1e-3;    // minimal value for first moment to avoid div by zero error
    for( unsigned j = 0; j < cellMids.size(); ++j ) {

        if( norm( cellMids[j] ) < 0.4 ) {
            kinetic_density = std::max( 0.2 * std::cos( norm( cellMids[j] ) * 4 ) * std::cos( norm( cellMids[j] ) * 4.0 ), epsilon );
        }
        else {
            kinetic_density = epsilon;
        }
        psi[j] = Vector( _settings->GetNQuadPoints(), kinetic_density );
    }
    return psi;
}

// ---- LineSource_PN ----

MeltingCube_Moment::MeltingCube_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}

MeltingCube_Moment::~MeltingCube_Moment() {}

VectorVector MeltingCube_Moment::SetupIC() {
    // Compute number of equations in the system

    // In case of PN, spherical basis is per default SPHERICAL_HARMONICS
    SphericalBase* tempBase  = SphericalBase::Create( _settings );
    unsigned ntotalEquations = tempBase->GetBasisSize();

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

    Vector uIC( ntotalEquations, 0 );

    if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
        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;
    }

    // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
    double epsilon = 1e-4;    // minimal value for first moment to avoid div by zero error
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double kinetic_density = 0.0;
        if( norm( cellMids[j] ) < 0.4 ) {
            kinetic_density = std::max( 0.2 * std::cos( norm( cellMids[j] ) * 4 ) * std::cos( norm( cellMids[j] ) * 4.0 ), epsilon );
        }
        else {
            kinetic_density = epsilon;
        }
        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            psi[j] = kinetic_density * uIC / uIC[0];    // Remember scaling
        }
        if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
            psi[j][0] = kinetic_density;
        }
    }
    delete tempBase;    // Only temporally needed
    return psi;
}

// ---- LineSource SN pseudo1D ----

MeltingCube_SN_1D::MeltingCube_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}

MeltingCube_SN_1D::~MeltingCube_SN_1D() {}

VectorVector MeltingCube_SN_1D::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    auto cellMids          = _mesh->GetCellMidPoints();
    double epsilon         = 1e-3;    // minimal value for first moment to avoid div by zero error
    double kinetic_density = 0.0;
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x = cellMids[j][0];
        if( x > -0.5 && x < 0 ) {
            kinetic_density = std::max( -x, epsilon );
        }
        else if( x > 0 && x < 0.5 ) {
            kinetic_density = std::max( -x, epsilon );
        }
        else {
            kinetic_density = epsilon;
        }
        psi[j] = kinetic_density;
    }
    return psi;
}

// ---- LineSource Moment pseudo1D ----

MeltingCube_Moment_1D::MeltingCube_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}

MeltingCube_Moment_1D::~MeltingCube_Moment_1D() {}

VectorVector MeltingCube_Moment_1D::SetupIC() {
    // double t     = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac impulse)
    double epsilon = 1e-3;    // minimal value for first moment to avoid div by zero error

    // In case of PN, spherical basis is per default SPHERICAL_HARMONICS
    if( _settings->GetSolverName() == PN_SOLVER || _settings->GetSolverName() == CSD_PN_SOLVER ) {
        // In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions
        SphericalHarmonics* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );
        unsigned ntotalEquations     = tempBase->GetBasisSize();
        delete tempBase;
        VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) );    // zero could lead to problems?
        VectorVector cellMids  = _mesh->GetCellMidPoints();
        double kinetic_density = 0.0;
        for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
            double x = cellMids[idx_cell][0];
            if( x > -0.4 && x < 0.4 ) {
                kinetic_density = std::max( 0.4 * std::cos( x * 4 ) * std::cos( x * 4.0 ), epsilon );
            }
            else {
                kinetic_density = epsilon;
            }
            initialSolution[idx_cell][0] = kinetic_density;
        }
        return initialSolution;
    }
    else {
        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 uIC( ntotalEquations, 0 );

        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            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;
        }
        // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
        for( unsigned j = 0; j < cellMids.size(); ++j ) {
            double x               = cellMids[j][0];
            double kinetic_density = 0.0;
            if( x > -0.4 && x < 0.4 ) {
                kinetic_density = std::max( 0.2 * std::cos( x * 4 ) * std::cos( x * 4.0 ), epsilon );
            }
            else {
                kinetic_density = epsilon;
            }
            if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
                initialSolution[j] = kinetic_density * uIC / uIC[0];    // Remember scaling
            }
            if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                initialSolution[j][0] = kinetic_density;
            }
        }
        delete tempBase;    // Only temporally needed
        return initialSolution;
    }
}

1	#include "problems/meltingcube.hpp"
2	#include "common/config.hpp"
3	#include "common/mesh.hpp"
4	#include "quadratures/quadraturebase.hpp"
5	#include "velocitybasis/sphericalbase.hpp"
6	#include "velocitybasis/sphericalharmonics.hpp"
7	#include <complex>
8
9	// ---- Linesource ----
10
NEW 11	MeltingCube::MeltingCube( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {	×
NEW 12	_sigmaS = settings->GetSigmaS();	×
13	}
14
15	MeltingCube::~MeltingCube() {}	×
16		×
17	VectorVector MeltingCube::GetScatteringXS( const Vector& energies ) {	×
18	return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) );
19	}	×
20		×
21	VectorVector MeltingCube::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) ); }
22
23	std::vector<VectorVector> MeltingCube::GetExternalSource( const Vector& /energies/ ) {	×
24	return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
25	}	×
26		×
27	// ---- MeltingCube_SN ----
28
29	MeltingCube_SN::MeltingCube_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}
30
31	MeltingCube_SN::~MeltingCube_SN() {}	×
32
33	VectorVector MeltingCube_SN::SetupIC() {	×
34	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );	×
35	auto cellMids = _mesh->GetCellMidPoints();	×
36	double kinetic_density = 0.0;
37	double epsilon = 1e-3; // minimal value for first moment to avoid div by zero error	×
38	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
NEW 39		×
40	if( norm( cellMids[j] ) < 0.4 ) {	×
41	kinetic_density = std::max( 0.2 * std::cos( norm( cellMids[j] ) * 4 ) * std::cos( norm( cellMids[j] ) * 4.0 ), epsilon );	×
42	}	×
43	else {
44	kinetic_density = epsilon;	×
45	}	×
46	psi[j] = Vector( _settings->GetNQuadPoints(), kinetic_density );
47	}
48	return psi;	×
49	}
50		×
51	// ---- LineSource_PN ----
52		×
53	MeltingCube_Moment::MeltingCube_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}
54
55	MeltingCube_Moment::~MeltingCube_Moment() {}
56
57	VectorVector MeltingCube_Moment::SetupIC() {	×
58	// Compute number of equations in the system
59		×
60	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
61	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
62	unsigned ntotalEquations = tempBase->GetBasisSize();
63		×
64	VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
65	VectorVector cellMids = _mesh->GetCellMidPoints();
66
67	Vector uIC( ntotalEquations, 0 );	×
68		×
69	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
70	QuadratureBase* quad = QuadratureBase::Create( _settings );	×
71	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
72	Vector w = quad->GetWeights();
73		×
74	double my, phi;
75	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
76		×
77	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
78	my = quadPointsSphere[idx_quad][0];	×
79	phi = quadPointsSphere[idx_quad][1];
80	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
81	}	×
82	// Integrate <1*m> to get factors for monomial basis in isotropic scattering
83	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
84	uIC += w[idx_quad] * moments[idx_quad];	×
85	}	×
86	delete quad;	×
87	}
88
89	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
90	double epsilon = 1e-4; // minimal value for first moment to avoid div by zero error	×
91	for( unsigned j = 0; j < cellMids.size(); ++j ) {
92	double kinetic_density = 0.0;	×
93	if( norm( cellMids[j] ) < 0.4 ) {
94	kinetic_density = std::max( 0.2 * std::cos( norm( cellMids[j] ) * 4 ) * std::cos( norm( cellMids[j] ) * 4.0 ), epsilon );
95	}
96	else {	×
97	kinetic_density = epsilon;	×
98	}	×
99	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
100	psi[j] = kinetic_density * uIC / uIC[0]; // Remember scaling	×
101	}
102	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
103	psi[j][0] = kinetic_density;	×
104	}
105	}	×
106	delete tempBase; // Only temporally needed	×
107	return psi;
108	}	×
109		×
110	// ---- LineSource SN pseudo1D ----
111
NEW 112	MeltingCube_SN_1D::MeltingCube_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}	×
113		×
114	MeltingCube_SN_1D::~MeltingCube_SN_1D() {}
115
116	VectorVector MeltingCube_SN_1D::SetupIC() {
117	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
118	auto cellMids = _mesh->GetCellMidPoints();	×
119	double epsilon = 1e-3; // minimal value for first moment to avoid div by zero error
120	double kinetic_density = 0.0;	×
121	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
122	double x = cellMids[j][0];	×
123	if( x > -0.5 && x < 0 ) {
124	kinetic_density = std::max( -x, epsilon );	×
125	}	×
126	else if( x > 0 && x < 0.5 ) {	×
127	kinetic_density = std::max( -x, epsilon );	×
128	}	×
129	else {	×
130	kinetic_density = epsilon;	×
131	}	×
132	psi[j] = kinetic_density;	×
133	}
134	return psi;	×
135	}	×
136
137	// ---- LineSource Moment pseudo1D ----
138		×
139	MeltingCube_Moment_1D::MeltingCube_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : MeltingCube( settings, mesh, quad ) {}
140		×
141	MeltingCube_Moment_1D::~MeltingCube_Moment_1D() {}
142		×
143	VectorVector MeltingCube_Moment_1D::SetupIC() {
144	// double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
145	double epsilon = 1e-3; // minimal value for first moment to avoid div by zero error
146
147	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
148	if( _settings->GetSolverName() == PN_SOLVER \|\| _settings->GetSolverName() == CSD_PN_SOLVER ) {
149	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions	×
150	SphericalHarmonics* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );	×
151	unsigned ntotalEquations = tempBase->GetBasisSize();	×
152	delete tempBase;
153	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?	×
154	VectorVector cellMids = _mesh->GetCellMidPoints();
155	double kinetic_density = 0.0;	×
156	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
157	double x = cellMids[idx_cell][0];
158	if( x > -0.4 && x < 0.4 ) {	×
159	kinetic_density = std::max( 0.4 * std::cos( x * 4 ) * std::cos( x * 4.0 ), epsilon );
160	}	×
161	else {	×
162	kinetic_density = epsilon;	×
163	}	×
164	initialSolution[idx_cell][0] = kinetic_density;	×
165	}	×
166	return initialSolution;	×
167	}	×
168	else {	×
169	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
170	unsigned ntotalEquations = tempBase->GetBasisSize();
171
172	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?	×
173	VectorVector cellMids = _mesh->GetCellMidPoints();
174	Vector uIC( ntotalEquations, 0 );	×
175
176	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
177	QuadratureBase* quad = QuadratureBase::Create( _settings );
178	VectorVector quadPointsSphere = quad->GetPointsSphere();
179	Vector w = quad->GetWeights();	×
180		×
181	double my, phi;
182	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
183		×
184	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
185	my = quadPointsSphere[idx_quad][0];
186	phi = quadPointsSphere[idx_quad][1];	×
187	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );	×
188	}	×
189	// Integrate <1*m> to get factors for monomial basis in isotropic scattering	×
190	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
191	uIC += w[idx_quad] * moments[idx_quad];
192	}	×
193	delete quad;
194	}	×
195	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
196	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
197	double x = cellMids[j][0];	×
198	double kinetic_density = 0.0;
199	if( x > -0.4 && x < 0.4 ) {
200	kinetic_density = std::max( 0.2 * std::cos( x * 4 ) * std::cos( x * 4.0 ), epsilon );	×
201	}	×
202	else {
203	kinetic_density = epsilon;	×
204	}
205	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
206	initialSolution[j] = kinetic_density * uIC / uIC[0]; // Remember scaling	×
207	}	×
208	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
209	initialSolution[j][0] = kinetic_density;	×
210	}	×
211	}
212	delete tempBase; // Only temporally needed
213	return initialSolution;	×
214	}
215	}	×

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