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

38.87

/src/problems/checkerboard.cpp

#include "problems/checkerboard.hpp"
#include "common/config.hpp"
#include "common/mesh.hpp"
#include "quadratures/quadraturebase.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"
#include "velocitybasis/sphericalbase.hpp"
#include "velocitybasis/sphericalharmonics.hpp"

// ---- Checkerboard Sn ----
// Constructor for Ckeckerboard case with Sn
Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {

    // Initialise crosssections to 1
    _sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
    _sigmaT = Vector( _mesh->GetNumCells(), 1.0 );

    // For absorption cells: set scattering XS to 0 and absorption to 10
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isAbsorption( cellMids[j] ) ) {
            _sigmaS[j] = 0.0;
            _sigmaT[j] = 10.0;
        }
    }
}

Checkerboard_SN::~Checkerboard_SN() {}

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

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

std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& /*energies*/ ) {
    VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isSource( cellMids[j] ) ) Q[j] = _settings->GetSourceMagnitude();    // isotropic source
    }
    return std::vector<VectorVector>( 1u, Q );
}

VectorVector Checkerboard_SN::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    return psi;
}

bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
    // Check whether pos is inside absorbing squares
    std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
    std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
    for( unsigned k = 0; k < lbounds.size(); ++k ) {
        for( unsigned l = 0; l < lbounds.size(); ++l ) {
            if( ( l + k ) % 2 == 1 || ( k == 2 && l == 2 ) || ( k == 2 && l == 4 ) ) continue;
            if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {
                return true;
            }
        }
    }
    return false;
}

bool Checkerboard_SN::isSource( const Vector& pos ) const {
    // Check whether pos is part of source region
    if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
        return true;
    else
        return false;
}

// ---- Checkerboard Moments ----

// Constructor for checkerboard case with Pn
Checkerboard_Moment::Checkerboard_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {

    // Initialise crosssections = 1 (scattering)
    _sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
    _sigmaT = Vector( _mesh->GetNumCells(), 1.0 );

    // for absorption regions change crosssections to all absorption
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isAbsorption( cellMids[j] ) ) {
            _sigmaS[j] = 0.0;
            _sigmaT[j] = 10.0;
        }
    }
}

Checkerboard_Moment::~Checkerboard_Moment() {}

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

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

std::vector<VectorVector> Checkerboard_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 ) {
        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 = _settings->GetSourceMagnitude();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isSource( cellMids[j] ) ) {
            if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
                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 Checkerboard_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 ) {
        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;
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
            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;
}

bool Checkerboard_Moment::isAbsorption( const Vector& pos ) const {
    // Check whether pos is in absorption region
    std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
    std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
    for( unsigned k = 0; k < lbounds.size(); ++k ) {
        for( unsigned l = 0; l < lbounds.size(); ++l ) {
            if( ( l + k ) % 2 == 1 || ( k == 2 && l == 2 ) || ( k == 2 && l == 4 ) ) continue;
            if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {
                return true;
            }
        }
    }
    return false;
}

bool Checkerboard_Moment::isSource( const Vector& pos ) const {
    // Check whether pos is in source region
    if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
        return true;
    else
        return false;
}

// ---- Checkerboard SN 1D ----
// Constructor for Ckeckerboard case with Sn
Checkerboard_SN_1D::Checkerboard_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {

    // Initialise crosssections to 1
    _sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
    _sigmaT = Vector( _mesh->GetNumCells(), 1.0 );

    // For absorption cells: set scattering XS to 0 and absorption to 10
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isAbsorption( cellMids[j] ) ) {
            _sigmaS[j] = 0.0;
            _sigmaT[j] = 10.0;
        }
    }
}

Checkerboard_SN_1D::~Checkerboard_SN_1D() {}

VectorVector Checkerboard_SN_1D::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _sigmaS ); }

VectorVector Checkerboard_SN_1D::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _sigmaT ); }

std::vector<VectorVector> Checkerboard_SN_1D::GetExternalSource( const Vector& /*energies*/ ) {
    VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isSource( cellMids[j] ) ) Q[j] = _settings->GetSourceMagnitude();    // isotropic source
    }
    return std::vector<VectorVector>( 1u, Q );
}

VectorVector Checkerboard_SN_1D::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    return psi;
}

bool Checkerboard_SN_1D::isAbsorption( const Vector& pos ) const {
    // Check whether pos is inside absorbing squares
    // domain from 0 to 7, absorption block is between 1 and 2
    if( ( pos[0] >= 1 && pos[0] <= 2 ) || ( pos[0] >= 6.5 && pos[0] <= 7 ) ) {
        return true;
    }
    return false;
}

bool Checkerboard_SN_1D::isSource( const Vector& pos ) const {
    // Check whether pos is part of source region
    if( pos[0] >= 3 && pos[0] <= 4 )
        return true;
    else
        return false;
}

// --- Moment version 1d ---

Checkerboard_Moment_1D::Checkerboard_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {

    // Initialise crosssections to 1
    _sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
    _sigmaT = Vector( _mesh->GetNumCells(), 1.0 );

    // For absorption cells: set scattering XS to 0 and absorption to 10
    auto cellMids = _mesh->GetCellMidPoints();
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        if( isAbsorption( cellMids[j] ) ) {
            _sigmaS[j] = 0.0;
            _sigmaT[j] = 10.0;
        }
    }
}

Checkerboard_Moment_1D::~Checkerboard_Moment_1D() {}

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

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

std::vector<VectorVector> Checkerboard_Moment_1D::GetExternalSource( const Vector& /*energies*/ ) {
    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 Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) );
        double kinetic_density = _settings->GetSourceMagnitude();
        VectorVector cellMids  = _mesh->GetCellMidPoints();
        for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
            if( isSource( cellMids[idx_cell] ) ) {
                Q[idx_cell][0] = kinetic_density;
            }
        }
        return std::vector<VectorVector>( 1u, Q );
    }
    else {
        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 ) {
            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 = _settings->GetSourceMagnitude();
        for( unsigned j = 0; j < cellMids.size(); ++j ) {
            if( isSource( cellMids[j] ) ) {
                if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
                    Q[j] = kinetic_density * uIC / uIC[0];    // Remember scaling
                }
                if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                    Q[j][0] = kinetic_density;
                }
            }
        }
        delete tempBase;    // Only temporally needed
        return std::vector<VectorVector>( 1u, Q );
    }
}

VectorVector Checkerboard_Moment_1D::SetupIC() {
    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
        SphericalBase* tempBase  = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );
        unsigned ntotalEquations = tempBase->GetBasisSize();
        delete tempBase;
        double epsilon = 1e-3;
        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 ) {
            initialSolution[idx_cell][0] = epsilon;
        }
        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 tempIC( 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++ ) {
                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;
        for( unsigned j = 0; j < cellMids.size(); ++j ) {
            if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
                initialSolution[j] = kinetic_density * tempIC / tempIC[0];    // Remember scaling
            }
            if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                initialSolution[j][0] = kinetic_density;
            }
        }
        delete tempBase;    // Only temporally needed
        return initialSolution;
    }
}

bool Checkerboard_Moment_1D::isAbsorption( const Vector& pos ) const {
    // Check whether pos is inside absorbing squares
    // domain from 0 to 7, absorption block is between 1 and 2
    if( ( pos[0] >= 1 && pos[0] <= 2 ) || ( pos[0] >= 6.5 && pos[0] <= 7 ) ) {
        return true;
    }
    return false;
}

bool Checkerboard_Moment_1D::isSource( const Vector& pos ) const {
    // Check whether pos is part of source region
    if( pos[0] >= 3 && pos[0] <= 4 )
        return true;
    else
        return false;
}

1	#include "problems/checkerboard.hpp"
2	#include "common/config.hpp"
3	#include "common/mesh.hpp"
4	#include "quadratures/quadraturebase.hpp"
5	#include "toolboxes/textprocessingtoolbox.hpp"
6	#include "velocitybasis/sphericalbase.hpp"
7	#include "velocitybasis/sphericalharmonics.hpp"
8
9	// ---- Checkerboard Sn ----
10	// Constructor for Ckeckerboard case with Sn
11	Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {	2✔
12
13	// Initialise crosssections to 1
14	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );	2✔
15	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );	2✔
16
17	// For absorption cells: set scattering XS to 0 and absorption to 10
18	auto cellMids = _mesh->GetCellMidPoints();	4✔
19	for( unsigned j = 0; j < cellMids.size(); ++j ) {	2,114✔
20	if( isAbsorption( cellMids[j] ) ) {	2,112✔
21	_sigmaS[j] = 0.0;	572✔
22	_sigmaT[j] = 10.0;	572✔
23	}
24	}
25	}	2✔
26
27	Checkerboard_SN::~Checkerboard_SN() {}	4✔
28		2✔
29	VectorVector Checkerboard_SN::GetScatteringXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaS ); }	2✔
30
31	VectorVector Checkerboard_SN::GetTotalXS( const Vector& /energies / ) { return VectorVector( 1u, _sigmaT ); }	2✔
32
33	std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& /energies/ ) {	2✔
34	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
35	auto cellMids = _mesh->GetCellMidPoints();	2✔
36	for( unsigned j = 0; j < cellMids.size(); ++j ) {	6✔
37	if( isSource( cellMids[j] ) ) Q[j] = _settings->GetSourceMagnitude(); // isotropic source	2✔
38	}	2,114✔
39	return std::vector<VectorVector>( 1u, Q );	2,112✔
40	}
41		4✔
42	VectorVector Checkerboard_SN::SetupIC() {
43	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
44	return psi;	2✔
45	}	4✔
46		2✔
47	bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
48	// Check whether pos is inside absorbing squares
49	std::vector<double> lbounds{ 1, 2, 3, 4, 5 };	2,112✔
50	std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
51	for( unsigned k = 0; k < lbounds.size(); ++k ) {	4,224✔
52	for( unsigned l = 0; l < lbounds.size(); ++l ) {	4,224✔
53	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	10,956✔
54	if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {	54,676✔
55	return true;	45,832✔
56	}	20,372✔
57	}	572✔
58	}
59	return false;
60	}
61		1,540✔
62	bool Checkerboard_SN::isSource( const Vector& pos ) const {
63	// Check whether pos is part of source region
64	if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )	2,112✔
65	return true;
66	else	2,112✔
67	return false;	52✔
68	}
69		2,060✔
70	// ---- Checkerboard Moments ----
71
72	// Constructor for checkerboard case with Pn
73	Checkerboard_Moment::Checkerboard_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
74
75	// Initialise crosssections = 1 (scattering)	4✔
76	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
77	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );
78		4✔
79	// for absorption regions change crosssections to all absorption	4✔
80	auto cellMids = _mesh->GetCellMidPoints();
81	for( unsigned j = 0; j < cellMids.size(); ++j ) {
82	if( isAbsorption( cellMids[j] ) ) {	8✔
83	_sigmaS[j] = 0.0;	4,228✔
84	_sigmaT[j] = 10.0;	4,224✔
85	}	1,144✔
86	}	1,144✔
87	}
88
89	Checkerboard_Moment::~Checkerboard_Moment() {}	4✔
90
91	VectorVector Checkerboard_Moment::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }	8✔
92		4✔
93	VectorVector Checkerboard_Moment::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaT ); }	4✔
94
95	std::vector<VectorVector> Checkerboard_Moment::GetExternalSource( const Vector& /energies/ ) {	4✔
96	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS
97		4✔
98	double integrationFactor = ( 4 * M_PI );
99	if( _settings->GetDim() == 2 ) {	4✔
100	integrationFactor = M_PI;
101	}
102	SphericalBase* tempBase = SphericalBase::Create( _settings );	4✔
103	unsigned ntotalEquations = tempBase->GetBasisSize();	4✔
104		×
105	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
106	VectorVector cellMids = _mesh->GetCellMidPoints();	4✔
107		4✔
108	Vector uIC( ntotalEquations, 0 );
109		12✔
110	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	8✔
111	QuadratureBase* quad = QuadratureBase::Create( _settings );
112	VectorVector quadPointsSphere = quad->GetPointsSphere();	4✔
113	Vector w = quad->GetWeights();
114		4✔
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 = _settings->GetSourceMagnitude();	×
130	for( unsigned j = 0; j < cellMids.size(); ++j ) {
131	if( isSource( cellMids[j] ) ) {	×
132	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
133	Q[j] = kinetic_density * uIC / uIC[0] / integrationFactor; // Remember scaling	4✔
134	}	4,228✔
135	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	4,224✔
136	Q[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	104✔
137	}	×
138	}
139	}	104✔
140	delete tempBase; // Only temporally needed	104✔
141	return std::vector<VectorVector>( 1u, Q );
142	}
143
144	VectorVector Checkerboard_Moment::SetupIC() {	4✔
145	double integrationFactor = ( 4 * M_PI );	8✔
146	if( _settings->GetDim() == 2 ) {
147	integrationFactor = M_PI;
148	}	4✔
149	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS	4✔
150	SphericalBase* tempBase = SphericalBase::Create( _settings );	4✔
151	unsigned ntotalEquations = tempBase->GetBasisSize();	×
152
153	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
154	VectorVector cellMids = _mesh->GetCellMidPoints();	4✔
155		4✔
156	Vector tempIC( ntotalEquations, 0 );
157		8✔
158	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	8✔
159	QuadratureBase* quad = QuadratureBase::Create( _settings );
160	VectorVector quadPointsSphere = quad->GetPointsSphere();	8✔
161	Vector w = quad->GetWeights();
162		4✔
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	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
180	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
181	initialSolution[j] = kinetic_density * tempIC / tempIC[0] / integrationFactor; // Remember scaling
182	}	4✔
183	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	4,228✔
184	initialSolution[j][0] = kinetic_density / integrationFactor; // first bassis function is 1/ ( 4 * M_PI )	4,224✔
185	}	×
186	}
187	delete tempBase; // Only temporally needed	4,224✔
188	return initialSolution;	4,224✔
189	}
190
191	bool Checkerboard_Moment::isAbsorption( const Vector& pos ) const {	4✔
192	// Check whether pos is in absorption region	8✔
193	std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
194	std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
195	for( unsigned k = 0; k < lbounds.size(); ++k ) {	4,224✔
196	for( unsigned l = 0; l < lbounds.size(); ++l ) {
197	if( ( l + k ) % 2 == 1 \|\| ( k == 2 && l == 2 ) \|\| ( k == 2 && l == 4 ) ) continue;	8,448✔
198	if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {	8,448✔
199	return true;	21,912✔
200	}	109,352✔
201	}	91,664✔
202	}	40,744✔
203	return false;	1,144✔
204	}
205
206	bool Checkerboard_Moment::isSource( const Vector& pos ) const {
207	// Check whether pos is in source region	3,080✔
208	if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
209	return true;
210	else	4,224✔
211	return false;
212	}	4,224✔
213		104✔
214	// ---- Checkerboard SN 1D ----
215	// Constructor for Ckeckerboard case with Sn	4,120✔
216	Checkerboard_SN_1D::Checkerboard_SN_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
217
218	// Initialise crosssections to 1
219	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
NEW 220	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );	×
221
222	// For absorption cells: set scattering XS to 0 and absorption to 10
223	auto cellMids = _mesh->GetCellMidPoints();	×
224	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
225	if( isAbsorption( cellMids[j] ) ) {
226	_sigmaS[j] = 0.0;
NEW 227	_sigmaT[j] = 10.0;	×
228	}	×
229	}	×
230	}	×
231		×
232	Checkerboard_SN_1D::~Checkerboard_SN_1D() {}
233
234	VectorVector Checkerboard_SN_1D::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _sigmaS ); }
235
NEW 236	VectorVector Checkerboard_SN_1D::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _sigmaT ); }	×
237		×
238	std::vector<VectorVector> Checkerboard_SN_1D::GetExternalSource( const Vector& /energies/ ) {	×
239	VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
240	auto cellMids = _mesh->GetCellMidPoints();	×
241	for( unsigned j = 0; j < cellMids.size(); ++j ) {
242	if( isSource( cellMids[j] ) ) Q[j] = _settings->GetSourceMagnitude(); // isotropic source	×
243	}
244	return std::vector<VectorVector>( 1u, Q );	×
245	}	×
246		×
247	VectorVector Checkerboard_SN_1D::SetupIC() {	×
248	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );	×
249	return psi;
250	}	×
251
252	bool Checkerboard_SN_1D::isAbsorption( const Vector& pos ) const {
253	// Check whether pos is inside absorbing squares	×
254	// domain from 0 to 7, absorption block is between 1 and 2	×
255	if( ( pos[0] >= 1 && pos[0] <= 2 ) \|\| ( pos[0] >= 6.5 && pos[0] <= 7 ) ) {	×
256	return true;
257	}
258	return false;	×
259	}
260
261	bool Checkerboard_SN_1D::isSource( const Vector& pos ) const {	×
262	// Check whether pos is part of source region	×
263	if( pos[0] >= 3 && pos[0] <= 4 )
264	return true;	×
265	else
266	return false;
267	}	×
268
269	// --- Moment version 1d ---	×
270		×
271	Checkerboard_Moment_1D::Checkerboard_Moment_1D( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {
272		×
273	// Initialise crosssections to 1
274	_sigmaS = Vector( _mesh->GetNumCells(), 1.0 );
275	_sigmaT = Vector( _mesh->GetNumCells(), 1.0 );
276
277	// For absorption cells: set scattering XS to 0 and absorption to 10	×
278	auto cellMids = _mesh->GetCellMidPoints();
279	for( unsigned j = 0; j < cellMids.size(); ++j ) {
280	if( isAbsorption( cellMids[j] ) ) {	×
NEW 281	_sigmaS[j] = 0.0;	×
282	_sigmaT[j] = 10.0;
283	}
284	}	×
285	}	×
286		×
287	Checkerboard_Moment_1D::~Checkerboard_Moment_1D() {}	×
288		×
289	VectorVector Checkerboard_Moment_1D::GetScatteringXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaS ); }
290
291	VectorVector Checkerboard_Moment_1D::GetTotalXS( const Vector& /energies/ ) { return VectorVector( 1u, _sigmaT ); }
292
293	std::vector<VectorVector> Checkerboard_Moment_1D::GetExternalSource( const Vector& /energies/ ) {	×
294	if( _settings->GetSolverName() == PN_SOLVER \|\| _settings->GetSolverName() == CSD_PN_SOLVER ) {	×
295	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions	×
296	SphericalHarmonics* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );
297	unsigned ntotalEquations = tempBase->GetBasisSize();	×
298	delete tempBase;
299	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) );	×
300	double kinetic_density = _settings->GetSourceMagnitude();
301	VectorVector cellMids = _mesh->GetCellMidPoints();	×
302	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {	×
303	if( isSource( cellMids[idx_cell] ) ) {
304	Q[idx_cell][0] = kinetic_density;	×
305	}	×
306	}	×
307	return std::vector<VectorVector>( 1u, Q );	×
308	}	×
309	else {	×
310	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
311	unsigned ntotalEquations = tempBase->GetBasisSize();	×
312		×
313	VectorVector Q( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
314	VectorVector cellMids = _mesh->GetCellMidPoints();
315	Vector uIC( ntotalEquations, 0 );	×
316	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
317	QuadratureBase* quad = QuadratureBase::Create( _settings );
318	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
319	Vector w = quad->GetWeights();	×
320
321	double my, phi;	×
322	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
323		×
324	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
325	my = quadPointsSphere[idx_quad][0];	×
326	phi = quadPointsSphere[idx_quad][1];	×
327	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );	×
328	}
329	// Integrate <1*m> to get factors for monomial basis in isotropic scattering
330	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
331	uIC += w[idx_quad] * moments[idx_quad];
332	}	×
333	delete quad;	×
334	}	×
335	double kinetic_density = _settings->GetSourceMagnitude();	×
336	for( unsigned j = 0; j < cellMids.size(); ++j ) {
337	if( isSource( cellMids[j] ) ) {
338	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
339	Q[j] = kinetic_density * uIC / uIC[0]; // Remember scaling	×
340	}
341	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
342	Q[j][0] = kinetic_density;
343	}	×
344	}	×
345	}	×
346	delete tempBase; // Only temporally needed	×
347	return std::vector<VectorVector>( 1u, Q );	×
348	}
349	}	×
350		×
351	VectorVector Checkerboard_Moment_1D::SetupIC() {
352	if( _settings->GetSolverName() == PN_SOLVER \|\| _settings->GetSolverName() == CSD_PN_SOLVER ) {
353	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions
354	SphericalBase* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );	×
355	unsigned ntotalEquations = tempBase->GetBasisSize();	×
356	delete tempBase;
357	double epsilon = 1e-3;
358	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0.0 ) ); // zero could lead to problems?
359	VectorVector cellMids = _mesh->GetCellMidPoints();	×
360	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {	×
361	initialSolution[idx_cell][0] = epsilon;
362	}	×
363	return initialSolution;	×
364	}	×
365	else {	×
366	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
367	unsigned ntotalEquations = tempBase->GetBasisSize();	×
368		×
369	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?	×
370	VectorVector cellMids = _mesh->GetCellMidPoints();
371	Vector tempIC( ntotalEquations, 0 );	×
372	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
373	QuadratureBase* quad = QuadratureBase::Create( _settings );
374	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
375	Vector w = quad->GetWeights();	×
376
377	double my, phi;	×
378	VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );	×
379		×
380	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
381	my = quadPointsSphere[idx_quad][0];	×
382	phi = quadPointsSphere[idx_quad][1];	×
383	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );	×
384	}
385	// Integrate <1*m> to get factors for monomial basis in isotropic scattering
386	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {	×
387	tempIC += w[idx_quad] * moments[idx_quad];
388	}	×
389	delete quad;	×
390	}	×
391	// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.	×
392	double kinetic_density = 1e-4;
393	for( unsigned j = 0; j < cellMids.size(); ++j ) {
394	if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {	×
395	initialSolution[j] = kinetic_density * tempIC / tempIC[0]; // Remember scaling	×
396	}
397	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
398	initialSolution[j][0] = kinetic_density;
399	}
400	}	×
401	delete tempBase; // Only temporally needed	×
402	return initialSolution;	×
403	}	×
404	}
405		×
406	bool Checkerboard_Moment_1D::isAbsorption( const Vector& pos ) const {	×
407	// Check whether pos is inside absorbing squares
408	// domain from 0 to 7, absorption block is between 1 and 2
409	if( ( pos[0] >= 1 && pos[0] <= 2 ) \|\| ( pos[0] >= 6.5 && pos[0] <= 7 ) ) {	×
410	return true;	×
411	}
412	return false;
413	}
414		×
415	bool Checkerboard_Moment_1D::isSource( const Vector& pos ) const {
416	// Check whether pos is part of source region
417	if( pos[0] >= 3 && pos[0] <= 4 )	×
418	return true;	×
419	else
420	return false;	×
421	}

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