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

30.36

/src/problems/linesource.cpp

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

// ---- Linesource ----

LineSource::LineSource( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
    _sigmaS = Vector( _mesh->GetNumCells(), _settings->GetSigmaS() );
    _sigmaT = Vector( _mesh->GetNumCells(), _settings->GetSigmaS() );
}

LineSource::~LineSource() {}

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

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

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

double LineSource::GetAnalyticalSolution( double x, double y, double t, double /*sigma_s*/ ) {

    double solution = 0.0;
    double R        = sqrt( x * x + y * y );

    if( t > 0 ) {
        if( _sigmaS == 0.0 ) {
            if( ( t - R ) > 0 ) {
                solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
            }
        }
        else if( _sigmaS == 1.0 ) {
            double gamma = R / t;

            if( ( 1 - gamma ) > 0 ) {
                solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) ) + 2 * t * HelperIntRho_ptc( R, t );
            }
        }
        else {
            solution = 0.0;
        }
    }

    return M_PI * solution;    // Scaling of soolution ( 4 * M_PI ) *
}

double LineSource::HelperIntRho_ptc( double R, double t ) {

    int numsteps    = 100;
    double integral = 0;
    double gamma    = R / t;
    double omega    = 0;
    // integral is from 0 to  sqrt( 1 - gamma * gamma )
    double stepsize = sqrt( 1 - gamma * gamma ) / (double)numsteps;

    for( int i = 0; i < numsteps; i++ ) {
        omega = i * stepsize + 0.5 * stepsize;

        integral += stepsize * HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );
    }
    return integral;
}

double LineSource::HelperRho_ptc( double R, double t ) {
    double result = HelperRho_ptc1( R, t ) + HelperRho_ptc2( R, t );
    return result;
}

double LineSource::HelperRho_ptc1( double R, double t ) {
    double gamma  = R / t;
    double result = exp( -t ) / ( 4.0 * M_PI * R * t ) * log( ( 1 + gamma ) / ( 1 - gamma ) );
    return result;
}

double LineSource::HelperRho_ptc2( double R, double t ) {
    double gamma  = R / t;
    double result = 0;
    if( 1 - gamma > 0 ) {
        // Compute the integralpart with midpoint rule
        result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma ) * HelperIntRho_ptc2( t, gamma );
    }
    return result;
}

double LineSource::HelperIntRho_ptc2( double t, double gamma ) {
    double u        = 0;
    double integral = 0;
    std::complex<double> beta( 0, 0 );
    double q = 0;
    std::complex<double> complexPart( 0, 0 );
    std::complex<double> com_one( 0, 1 );

    // compute the integral using the midpoint rule
    // integral is from 0 to pi

    int numsteps = 100;

    double stepsize = M_PI / (double)numsteps;

    q = ( 1.0 + gamma ) / ( 1.0 - gamma );

    for( int i = 0; i < numsteps; i++ ) {
        u = i * stepsize + 0.5 * stepsize;

        // function evaluation
        beta        = ( log( q ) + com_one * u ) / ( gamma + com_one * tan( 0.5 * u ) );
        complexPart = ( gamma + com_one * tan( 0.5 * u ) ) * beta * beta * beta * exp( 0.5 * t * ( 1 - gamma * gamma ) * beta );
        integral += stepsize * ( 1 / cos( 0.5 * u ) ) * ( 1 / cos( 0.5 * u ) ) * complexPart.real();
    }
    return integral;
}

// ---- LineSource_SN ----

LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource( settings, mesh, quad ) {}

LineSource_SN::~LineSource_SN() {}

VectorVector LineSource_SN::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    auto cellMids = _mesh->GetCellMidPoints();
    double t      = 3.2e-4;    // pseudo time for gaussian smoothing
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x = cellMids[j][0];
        double y = cellMids[j][1];
        psi[j]   = Vector( _settings->GetNQuadPoints(), 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) ) / ( 4 * M_PI );
    }
    return psi;
}

// ---- LineSource_PN ----

LineSource_Moment::LineSource_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource( settings, mesh, quad ) {}

LineSource_Moment::~LineSource_Moment() {}

std::vector<VectorVector> LineSource_Moment::GetExternalSource( const Vector& /*energies*/ ) {
    SphericalBase* tempBase  = SphericalBase::Create( _settings );
    unsigned ntotalEquations = tempBase->GetBasisSize();
    delete tempBase;    // Only temporally needed

    return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ) );
}

VectorVector LineSource_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 initial_sol( _mesh->GetNumCells(), Vector( ntotalEquations, 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;
    }

    // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
    double t       = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac impulse)
    double epsilon = 1e-4;      // minimal value for first moment to avoid div by zero error

    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x = cellMids[j][0];
        double y = cellMids[j][1];    // (x- 0.5) * (x- 0.5)

        double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ), epsilon );

        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
            initial_sol[j] = kinetic_density * uIC / uIC[0];    // Remember scaling
        }
        if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
            initial_sol[j][0] = kinetic_density;
        }
    }
    delete tempBase;    // Only temporally needed
    return initial_sol;
}

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

LineSource_SN_1D::LineSource_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource_SN( settings, mesh, quad ) {}

VectorVector LineSource_SN_1D::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    auto cellMids  = _mesh->GetCellMidPoints();
    double t       = 3.2e-4;    // pseudo time for gaussian smoothing
    double epsilon = 1e-3;      // minimal value for first moment to avoid div by zero error

    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x               = cellMids[j][0];
        double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), epsilon );
        psi[j]                 = kinetic_density;
    }
    return psi;
}

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

LineSource_Moment_1D::LineSource_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource_Moment( settings, mesh, quad ) {}

VectorVector LineSource_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();

        for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
            double x                     = cellMids[idx_cell][0];
            double kinetic_density       = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), 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 || _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;
        }
        // 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 = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), epsilon );
            if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS || _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
                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/linesource.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
11	LineSource::LineSource( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {	12✔
12	_sigmaS = Vector( _mesh->GetNumCells(), _settings->GetSigmaS() );	12✔
13	_sigmaT = Vector( _mesh->GetNumCells(), _settings->GetSigmaS() );	12✔
14	}	12✔
15
16	LineSource::~LineSource() {}	10✔
17		×
18	VectorVector LineSource::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }	10✔
19
20	VectorVector LineSource::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }	12✔
21
22	std::vector<VectorVector> LineSource::GetExternalSource( const Vector& /energies/ ) {	12✔
23	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
24	return std::vector<VectorVector>( 1u, Q );	4✔
25	}	8✔
26		8✔
27	double LineSource::GetAnalyticalSolution( double x, double y, double t, double /sigma_s/ ) {
28
29	double solution = 0.0;	×
30	double R = sqrt( x * x + y * y );
31		×
32	if( t > 0 ) {	×
33	if( _sigmaS == 0.0 ) {
34	if( ( t - R ) > 0 ) {	×
35	solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );	×
36	}	×
37	}	×
38	else if( _sigmaS == 1.0 ) {
39	double gamma = R / t;
40		×
41	if( ( 1 - gamma ) > 0 ) {	×
42	solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) ) + 2 * t * HelperIntRho_ptc( R, t );
43	}	×
44	}	×
45	else {
46	solution = 0.0;
47	}
48	}	×
49
50	return M_PI * solution; // Scaling of soolution ( 4 * M_PI ) *
51	}
52		×
53	double LineSource::HelperIntRho_ptc( double R, double t ) {
54
55	int numsteps = 100;	×
56	double integral = 0;
57	double gamma = R / t;	×
58	double omega = 0;	×
59	// integral is from 0 to sqrt( 1 - gamma * gamma )	×
60	double stepsize = sqrt( 1 - gamma * gamma ) / (double)numsteps;	×
61
62	for( int i = 0; i < numsteps; i++ ) {	×
63	omega = i * stepsize + 0.5 * stepsize;
64		×
65	integral += stepsize * HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );	×
66	}
67	return integral;	×
68	}
69		×
70	double LineSource::HelperRho_ptc( double R, double t ) {
71	double result = HelperRho_ptc1( R, t ) + HelperRho_ptc2( R, t );
72	return result;	×
73	}	×
74		×
75	double LineSource::HelperRho_ptc1( double R, double t ) {
76	double gamma = R / t;
77	double result = exp( -t ) / ( 4.0 * M_PI * R * t ) * log( ( 1 + gamma ) / ( 1 - gamma ) );	×
78	return result;	×
79	}	×
80		×
81	double LineSource::HelperRho_ptc2( double R, double t ) {
82	double gamma = R / t;
83	double result = 0;	×
84	if( 1 - gamma > 0 ) {	×
85	// Compute the integralpart with midpoint rule	×
86	result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma ) * HelperIntRho_ptc2( t, gamma );	×
87	}
88	return result;	×
89	}
90		×
91	double LineSource::HelperIntRho_ptc2( double t, double gamma ) {
92	double u = 0;
93	double integral = 0;	×
94	std::complex<double> beta( 0, 0 );	×
95	double q = 0;	×
96	std::complex<double> complexPart( 0, 0 );	×
97	std::complex<double> com_one( 0, 1 );	×
98		×
99	// compute the integral using the midpoint rule	×
100	// integral is from 0 to pi
101
102	int numsteps = 100;
103
104	double stepsize = M_PI / (double)numsteps;	×
105
106	q = ( 1.0 + gamma ) / ( 1.0 - gamma );	×
107
108	for( int i = 0; i < numsteps; i++ ) {	×
109	u = i * stepsize + 0.5 * stepsize;
110		×
111	// function evaluation	×
112	beta = ( log( q ) + com_one * u ) / ( gamma + com_one * tan( 0.5 * u ) );
113	complexPart = ( gamma + com_one * tan( 0.5 * u ) ) * beta * beta * beta * exp( 0.5 * t * ( 1 - gamma * gamma ) * beta );
114	integral += stepsize * ( 1 / cos( 0.5 * u ) ) * ( 1 / cos( 0.5 * u ) ) * complexPart.real();	×
115	}	×
116	return integral;	×
117	}
118		×
119	// ---- LineSource_SN ----
120
121	LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource( settings, mesh, quad ) {}
122
123	LineSource_SN::~LineSource_SN() {}	4✔
124
125	VectorVector LineSource_SN::SetupIC() {	4✔
126	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );	2✔
127	auto cellMids = _mesh->GetCellMidPoints();	2✔
128	double t = 3.2e-4; // pseudo time for gaussian smoothing
129	for( unsigned j = 0; j < cellMids.size(); ++j ) {	4✔
130	double x = cellMids[j][0];	8✔
131	double y = cellMids[j][1];	8✔
132	psi[j] = Vector( _settings->GetNQuadPoints(), 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) ) / ( 4 * M_PI );	4✔
133	}	388✔
134	return psi;	384✔
135	}	384✔
136		384✔
137	// ---- LineSource_PN ----
138		8✔
139	LineSource_Moment::LineSource_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource( settings, mesh, quad ) {}
140
141	LineSource_Moment::~LineSource_Moment() {}
142
143	std::vector<VectorVector> LineSource_Moment::GetExternalSource( const Vector& /energies/ ) {	8✔
144	SphericalBase* tempBase = SphericalBase::Create( _settings );
145	unsigned ntotalEquations = tempBase->GetBasisSize();	16✔
146	delete tempBase; // Only temporally needed	8✔
147		8✔
148	return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ) );
149	}	8✔
150		8✔
151	VectorVector LineSource_Moment::SetupIC() {	8✔
152	// Compute number of equations in the system	8✔
153
154	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	16✔
155	SphericalBase* tempBase = SphericalBase::Create( _settings );
156	unsigned ntotalEquations = tempBase->GetBasisSize();
157		8✔
158	VectorVector initial_sol( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
159	VectorVector cellMids = _mesh->GetCellMidPoints();
160
161	Vector uIC( ntotalEquations, 0 );	8✔
162		8✔
163	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {
164	QuadratureBase* quad = QuadratureBase::Create( _settings );	16✔
165	VectorVector quadPointsSphere = quad->GetPointsSphere();	16✔
166	Vector w = quad->GetWeights();
167		16✔
168	double my, phi;
169	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	8✔
170		×
171	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
172	my = quadPointsSphere[idx_quad][0];	×
173	phi = quadPointsSphere[idx_quad][1];
174	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
175	}	×
176	// Integrate <1*m> to get factors for monomial basis in isotropic scattering
177	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
178	uIC += w[idx_quad] * moments[idx_quad];	×
179	}	×
180	delete quad;	×
181	}
182
183	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
184	double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)	×
185	double epsilon = 1e-4; // minimal value for first moment to avoid div by zero error
186		×
187	for( unsigned j = 0; j < cellMids.size(); ++j ) {
188	double x = cellMids[j][0];
189	double y = cellMids[j][1]; // (x- 0.5) * (x- 0.5)
190		8✔
191	double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ), epsilon );	8✔
192
193	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	776✔
194	initial_sol[j] = kinetic_density * uIC / uIC[0]; // Remember scaling	768✔
195	}	768✔
196	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
197	initial_sol[j][0] = kinetic_density;	768✔
198	}
199	}	768✔
200	delete tempBase; // Only temporally needed	×
201	return initial_sol;
202	}	768✔
203		768✔
204	// ---- LineSource SN pseudo1D ----
205
206	LineSource_SN_1D::LineSource_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource_SN( settings, mesh, quad ) {}	8✔
207		16✔
208	VectorVector LineSource_SN_1D::SetupIC() {
209	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
210	auto cellMids = _mesh->GetCellMidPoints();
211	double t = 3.2e-4; // pseudo time for gaussian smoothing
212	double epsilon = 1e-3; // minimal value for first moment to avoid div by zero error	×
213
214	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
215	double x = cellMids[j][0];	×
216	double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), epsilon );	×
217	psi[j] = kinetic_density;	×
218	}	×
219	return psi;
220	}	×
221		×
222	// ---- LineSource Moment pseudo1D ----	×
223		×
224	LineSource_Moment_1D::LineSource_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : LineSource_Moment( settings, mesh, quad ) {}
225		×
226	VectorVector LineSource_Moment_1D::SetupIC() {
227	double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
228	double epsilon = 1e-3; // minimal value for first moment to avoid div by zero error
229
230	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	×
231	if( _settings->GetSolverName() == PN_SOLVER \|\| _settings->GetSolverName() == CSD_PN_SOLVER ) {
232	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions	×
233	SphericalHarmonics* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );	×
234	unsigned ntotalEquations = tempBase->GetBasisSize();	×
235	delete tempBase;
236	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
237	VectorVector cellMids = _mesh->GetCellMidPoints();	×
238
239	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {	×
240	double x = cellMids[idx_cell][0];	×
241	double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), epsilon );	×
242	initialSolution[idx_cell][0] = kinetic_density;	×
243	}	×
244	return initialSolution;
245	}	×
246	else {	×
247	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
248	unsigned ntotalEquations = tempBase->GetBasisSize();	×
249
250	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?	×
251	VectorVector cellMids = _mesh->GetCellMidPoints();
252	Vector uIC( ntotalEquations, 0 );
253		×
NEW 254	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
255	QuadratureBase* quad = QuadratureBase::Create( _settings );
256	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
257	Vector w = quad->GetWeights();	×
258		×
259	double my, phi;
260	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
261		×
262	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
263	my = quadPointsSphere[idx_quad][0];	×
264	phi = quadPointsSphere[idx_quad][1];
265	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
266	}	×
267	// Integrate <1*m> to get factors for monomial basis in isotropic scattering
268	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
269	uIC += w[idx_quad] * moments[idx_quad];	×
270	}	×
271	delete quad;	×
272	}
273	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
274	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
275	double x = cellMids[j][0];	×
276	double kinetic_density = std::max( 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) ), epsilon );
NEW 277	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS \|\| _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS_ROTATED ) {	×
278	initialSolution[j] = kinetic_density * uIC / uIC[0]; // Remember scaling
279	}
280	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
281	initialSolution[j][0] = kinetic_density;	×
282	}	×
283	}	×
284	delete tempBase; // Only temporally needed	×
285	return initialSolution;
286	}	×
287	}	×

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