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

#include "problems/radiationctimage.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"

#include <fstream>

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

RadiationCTImage::~RadiationCTImage() {}

std::vector<VectorVector> RadiationCTImage::GetExternalSource( const Vector& energies ) {
    auto zeroVec = Vector( _settings->GetNQuadPoints(), 0.0 );
    auto uniform = std::vector<Vector>( _mesh->GetNumCells(), zeroVec );
    auto Q       = std::vector<VectorVector>( energies.size(), uniform );
    return Q;
}

VectorVector RadiationCTImage::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    VectorVector cellMids         = _mesh->GetCellMidPoints();
    double s                      = 0.1;
    double enterPositionX         = 2.5;    // 0.0;
    double enterPositionY         = 5.8;
    double meanDir                = M_PI / 2;
    double s_ang                  = 0.1;
    double epsilon                = 1e-3;
    QuadratureBase* quad          = QuadratureBase::Create( _settings );
    VectorVector quadPointsSphere = quad->GetPointsSphere();

    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x = cellMids[j][0] - enterPositionX;
        double y = cellMids[j][1] - enterPositionY;
        // anisotropic inflow that concentrates all particles on the last quadrature point
        //    for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
        //         double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]*quadPointsSphere[idx_quad][0])*cos(quadPointsSphere[idx_quad][1]));
        //           if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {
        //                 psi[j][idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) )  * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ), epsilon );
        //         }
        //         }
        // normal distribution also in angle

        for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
            if( quadPointsSphere[idx_quad][0] < 0 ) {
                double theta =
                    acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );
                double ang       = theta - meanDir;
                psi[j][idx_quad] = std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *
                                                 std::exp( -( ang * ang ) / ( 2 * s_ang ) ),
                                             epsilon );
            }
        }
    }

    delete quad;
    return psi;
}

std::vector<double> RadiationCTImage::GetDensity( const VectorVector& /*cellMidPoints*/ ) {
    std::string imageFile = _settings->GetCTFile();
    std::string meshFile  = _settings->GetMeshFile();
    ErrorMessages::Error( "Python API support is deprecated for KiT-RT. This function needs to be rewritten.\n Use a python script to first create "
                          "the mesh, then call KiT-RT from Python. This is the recommended workflow for stability and performance.",
                          CURRENT_FUNCTION );
    Matrix gsImage = Matrix( 0, 0 );    //    createSU2MeshFromImage( imageFile, meshFile ); DEprecated

    ErrorMessages::Error( "GetBounds function currently deprecetaded", CURRENT_FUNCTION );
    auto bounds        = _mesh->GetBounds();
    auto cellMidPoints = _mesh->GetCellMidPoints();

    double xMin = bounds[0].first;
    double xMax = bounds[0].second;
    double yMin = bounds[1].first;
    double yMax = bounds[1].second;

    unsigned m = gsImage.rows();
    unsigned n = gsImage.columns();

    Vector x( m ), y( n );
    for( unsigned i = 0; i < m; ++i ) {
        x[i] = xMin + static_cast<double>( i ) / static_cast<double>( m - 1 ) * ( xMax - xMin );
    }
    for( unsigned i = 0; i < n; ++i ) y[i] = yMin + static_cast<double>( i ) / static_cast<double>( n - 1 ) * ( yMax - yMin );

    Interpolation interp( x, y, gsImage );
    std::vector<double> result( _mesh->GetNumCells(), 0.0 );
    for( unsigned i = 0; i < _mesh->GetNumCells(); ++i ) {
        result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ) * 1.85,
                                0.05,
                                1.85 );    // Scale densities for CT to be between 0 (air) and 1.85 (bone)
    }
    return result;
}

VectorVector RadiationCTImage::GetScatteringXS( const Vector& /*energies*/ ) {
    // @TODO
    // Specified in subclasses
    return VectorVector( 1, Vector( 1, 0.0 ) );
}

VectorVector RadiationCTImage::GetTotalXS( const Vector& /*energies*/ ) {
    // @TODO
    // Specified in subclasses
    return VectorVector( 1, Vector( 1, 0.0 ) );
}

RadiationCTImage_Moment::RadiationCTImage_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : RadiationCTImage( settings, mesh, quad ) {}

RadiationCTImage_Moment::~RadiationCTImage_Moment() {}

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

    Vector zeroVec              = Vector( ntotalEquations, 0.0 );
    VectorVector uniform        = VectorVector( _mesh->GetNumCells(), zeroVec );
    std::vector<VectorVector> Q = std::vector<VectorVector>( energies.size(), uniform );
    return Q;
}

VectorVector RadiationCTImage_Moment::SetupIC() {
    if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
        // 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();

        double epsilon = 1e-3;

        VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );
        VectorVector cellMids = _mesh->GetCellMidPoints();

        QuadratureBase* quad          = new QGaussLegendreTensorized( _settings );
        VectorVector quadPointsSphere = quad->GetPointsSphere();
        Vector w                      = quad->GetWeights();
        Vector cellKineticDensity( quad->GetNq(), epsilon );

        // compute moment basis
        VectorVector moments = VectorVector( quad->GetNq(), Vector( ntotalEquations, 0.0 ) );
        double my, phi;    // quadpoints in spherical coordinates

        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 );
        }
        delete tempBase;
        double s              = 0.1;
        double enterPositionX = 2.5;    // 0.0;
        double enterPositionY = 5.8;
        double meanDir        = M_PI / 2;
        double s_ang          = 0.1;

        for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
            double x = cellMids[idx_cell][0] - enterPositionX;
            double y = cellMids[idx_cell][1] - enterPositionY;
            // anisotropic, forward-directed particle inflow
            // for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
            //     double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]*quadPointsSphere[idx_quad][0])*cos(quadPointsSphere[idx_quad][1]));
            //     if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {
            //              cellKineticDensity[idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) )  * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ),
            //              epsilon );
            //      }
            //  }

            // normal distribution also in angle

            for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
                if( quadPointsSphere[idx_quad][0] < 0 ) {
                    double theta =
                        acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );
                    double ang = theta - meanDir;
                    cellKineticDensity[idx_quad] =
                        std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *
                                      std::exp( -( ang * ang ) / ( 2 * s_ang ) ),
                                  epsilon );
                }
            }

            // Compute moments of this kinetic density
            for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
                initialSolution[idx_cell] += cellKineticDensity[idx_quad] * w[idx_quad] * moments[idx_quad];
            }
        }
        // TextProcessingToolbox::PrintVectorVector(initialSolution);
        // exit(0);
        delete quad;
        return initialSolution;
    }
    else {
        SphericalBase* tempBase  = SphericalBase::Create( _settings );
        unsigned ntotalEquations = tempBase->GetBasisSize();

        double epsilon = 1e-3;

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

        QuadratureBase* quad          = QuadratureBase::Create( _settings );
        VectorVector quadPointsSphere = quad->GetPointsSphere();
        Vector w                      = quad->GetWeights();

        // compute moment basis
        VectorVector moments = VectorVector( quad->GetNq(), Vector( ntotalEquations, 0.0 ) );
        double my, phi;    // quadpoints in spherical coordinates
        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 );
        }
        delete tempBase;

        double s              = 0.1;
        double enterPositionX = 2.5;    // 0.0;
        double enterPositionY = 5.8;
        double meanDir        = M_PI / 2;
        double s_ang          = 0.1;

        for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {
            double x = cellMids[idx_cell][0] - enterPositionX;
            double y = cellMids[idx_cell][1] - enterPositionY;
            // anisotropic, forward-directed particle inflow
            // for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
            //     double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]*quadPointsSphere[idx_quad][0])*cos(quadPointsSphere[idx_quad][1]));
            //     if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {
            //              cellKineticDensity[idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) )  * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ),
            //              epsilon );
            //      }
            //  }

            // normal distribution also in angle

            for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
                if( quadPointsSphere[idx_quad][0] < 0 ) {
                    double theta =
                        acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );
                    double ang = theta - meanDir;
                    cellKineticDensity[idx_quad] =
                        std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *
                                      std::exp( -( ang * ang ) / ( 2 * s_ang ) ),
                                  epsilon );
                }
            }

            // Compute moments of this kinetic density
            for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
                initialSolution[idx_cell] += cellKineticDensity[idx_quad] * w[idx_quad] * moments[idx_quad];
            }
        }
        delete quad;
        return initialSolution;
    }
}

1	#include "problems/radiationctimage.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
14	#include <fstream>
15
NEW 16	RadiationCTImage::RadiationCTImage( Config* settings, Mesh* mesh, QuadratureBase* quad ) : ProblemBase( settings, mesh, quad ) {}	×
17
18	RadiationCTImage::~RadiationCTImage() {}	×
19		×
20	std::vector<VectorVector> RadiationCTImage::GetExternalSource( const Vector& energies ) {	×
21	auto zeroVec = Vector( _settings->GetNQuadPoints(), 0.0 );
22	auto uniform = std::vector<Vector>( _mesh->GetNumCells(), zeroVec );	×
23	auto Q = std::vector<VectorVector>( energies.size(), uniform );	×
24	return Q;	×
25	}	×
26		×
27	VectorVector RadiationCTImage::SetupIC() {
28	VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
29	VectorVector cellMids = _mesh->GetCellMidPoints();	×
30	double s = 0.1;	×
31	double enterPositionX = 2.5; // 0.0;	×
32	double enterPositionY = 5.8;	×
33	double meanDir = M_PI / 2;	×
34	double s_ang = 0.1;	×
35	double epsilon = 1e-3;	×
36	QuadratureBase* quad = QuadratureBase::Create( _settings );	×
37	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
38		×
39	for( unsigned j = 0; j < cellMids.size(); ++j ) {	×
40	double x = cellMids[j][0] - enterPositionX;
41	double y = cellMids[j][1] - enterPositionY;	×
42	// anisotropic inflow that concentrates all particles on the last quadrature point	×
43	// for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {	×
44	// double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]quadPointsSphere[idx_quad][0])cos(quadPointsSphere[idx_quad][1]));
45	// if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {
46	// psi[j][idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ), epsilon );
47	// }
48	// }
49	// normal distribution also in angle
50
51	for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
52	if( quadPointsSphere[idx_quad][0] < 0 ) {
53	double theta =	×
54	acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );	×
55	double ang = theta - meanDir;
56	psi[j][idx_quad] = std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *	×
57	std::exp( -( ang * ang ) / ( 2 * s_ang ) ),	×
58	epsilon );	×
59	}	×
60	}	×
61	}
62
63	delete quad;
64	return psi;
65	}	×
66		×
67	std::vector<double> RadiationCTImage::GetDensity( const VectorVector& /cellMidPoints/ ) {
68	std::string imageFile = _settings->GetCTFile();
69	std::string meshFile = _settings->GetMeshFile();	×
NEW 70	ErrorMessages::Error( "Python API support is deprecated for KiT-RT. This function needs to be rewritten.\n Use a python script to first create "	×
NEW 71	"the mesh, then call KiT-RT from Python. This is the recommended workflow for stability and performance.",	×
NEW 72	CURRENT_FUNCTION );	×
73	Matrix gsImage = Matrix( 0, 0 ); // createSU2MeshFromImage( imageFile, meshFile ); DEprecated
74
NEW 75	ErrorMessages::Error( "GetBounds function currently deprecetaded", CURRENT_FUNCTION );	×
76	auto bounds = _mesh->GetBounds();
NEW 77	auto cellMidPoints = _mesh->GetCellMidPoints();	×
78		×
79	double xMin = bounds[0].first;	×
80	double xMax = bounds[0].second;
81	double yMin = bounds[1].first;	×
82	double yMax = bounds[1].second;	×
83		×
84	unsigned m = gsImage.rows();	×
85	unsigned n = gsImage.columns();
86		×
87	Vector x( m ), y( n );	×
88	for( unsigned i = 0; i < m; ++i ) {
89	x[i] = xMin + static_cast<double>( i ) / static_cast<double>( m - 1 ) * ( xMax - xMin );	×
90	}	×
91	for( unsigned i = 0; i < n; ++i ) y[i] = yMin + static_cast<double>( i ) / static_cast<double>( n - 1 ) * ( yMax - yMin );	×
92
93	Interpolation interp( x, y, gsImage );	×
94	std::vector<double> result( _mesh->GetNumCells(), 0.0 );
95	for( unsigned i = 0; i < _mesh->GetNumCells(); ++i ) {	×
96	result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ) * 1.85,	×
97	0.05,	×
98	1.85 ); // Scale densities for CT to be between 0 (air) and 1.85 (bone)	×
99	}	×
100	return result;	×
101	}
102		×
103	VectorVector RadiationCTImage::GetScatteringXS( const Vector& /energies/ ) {
104	// @TODO
105	// Specified in subclasses	×
106	return VectorVector( 1, Vector( 1, 0.0 ) );
107	}
108		×
109	VectorVector RadiationCTImage::GetTotalXS( const Vector& /energies/ ) {
110	// @TODO
111	// Specified in subclasses	×
112	return VectorVector( 1, Vector( 1, 0.0 ) );
113	}
114		×
115	RadiationCTImage_Moment::RadiationCTImage_Moment( Config* settings, Mesh* mesh, QuadratureBase* quad ) : RadiationCTImage( settings, mesh, quad ) {}
116
117	RadiationCTImage_Moment::~RadiationCTImage_Moment() {}	×
118
119	std::vector<VectorVector> RadiationCTImage_Moment::GetExternalSource( const Vector& energies ) {	×
120	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
121	unsigned ntotalEquations = tempBase->GetBasisSize();	×
122	delete tempBase;
123		×
124	Vector zeroVec = Vector( ntotalEquations, 0.0 );	×
125	VectorVector uniform = VectorVector( _mesh->GetNumCells(), zeroVec );	×
126	std::vector<VectorVector> Q = std::vector<VectorVector>( energies.size(), uniform );	×
127	return Q;
128	}	×
129		×
130	VectorVector RadiationCTImage_Moment::SetupIC() {	×
131	if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {	×
132	// In case of PN, spherical basis is per default SPHERICAL_HARMONICS in 3 velocity dimensions
133	SphericalBase* tempBase = new SphericalHarmonics( _settings->GetMaxMomentDegree(), 3 );
134	unsigned ntotalEquations = tempBase->GetBasisSize();	×
135		×
136	double epsilon = 1e-3;
137		×
138	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );	×
139	VectorVector cellMids = _mesh->GetCellMidPoints();
140		×
141	QuadratureBase* quad = new QGaussLegendreTensorized( _settings );
142	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
143	Vector w = quad->GetWeights();	×
144	Vector cellKineticDensity( quad->GetNq(), epsilon );
145		×
146	// compute moment basis	×
147	VectorVector moments = VectorVector( quad->GetNq(), Vector( ntotalEquations, 0.0 ) );	×
148	double my, phi; // quadpoints in spherical coordinates	×
149
150	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
151	my = quadPointsSphere[idx_quad][0];	×
152	phi = quadPointsSphere[idx_quad][1];
153	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
154	}	×
155	delete tempBase;	×
156	double s = 0.1;	×
157	double enterPositionX = 2.5; // 0.0;	×
158	double enterPositionY = 5.8;
159	double meanDir = M_PI / 2;	×
160	double s_ang = 0.1;	×
161		×
162	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {	×
163	double x = cellMids[idx_cell][0] - enterPositionX;	×
164	double y = cellMids[idx_cell][1] - enterPositionY;	×
165	// anisotropic, forward-directed particle inflow
166	// for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {	×
167	// double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]quadPointsSphere[idx_quad][0])cos(quadPointsSphere[idx_quad][1]));	×
168	// if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {	×
169	// cellKineticDensity[idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ),
170	// epsilon );
171	// }
172	// }
173
174	// normal distribution also in angle
175
176	for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
177	if( quadPointsSphere[idx_quad][0] < 0 ) {
178	double theta =
179	acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );
180	double ang = theta - meanDir;	×
181	cellKineticDensity[idx_quad] =	×
182	std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *
183	std::exp( -( ang * ang ) / ( 2 * s_ang ) ),	×
184	epsilon );	×
185	}	×
186	}	×
187		×
188	// Compute moments of this kinetic density	×
189	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
190	initialSolution[idx_cell] += cellKineticDensity[idx_quad] * w[idx_quad] * moments[idx_quad];
191	}
192	}
193	// TextProcessingToolbox::PrintVectorVector(initialSolution);	×
194	// exit(0);	×
195	delete quad;
196	return initialSolution;
197	}
198	else {
199	SphericalBase* tempBase = SphericalBase::Create( _settings );	×
200	unsigned ntotalEquations = tempBase->GetBasisSize();	×
201
202	double epsilon = 1e-3;
203		×
204	Vector cellKineticDensity( _settings->GetNQuadPoints(), epsilon );	×
205	VectorVector initialSolution( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
206	VectorVector cellMids = _mesh->GetCellMidPoints();	×
207
208	QuadratureBase* quad = QuadratureBase::Create( _settings );	×
209	VectorVector quadPointsSphere = quad->GetPointsSphere();	×
210	Vector w = quad->GetWeights();	×
211
212	// compute moment basis	×
213	VectorVector moments = VectorVector( quad->GetNq(), Vector( ntotalEquations, 0.0 ) );	×
214	double my, phi; // quadpoints in spherical coordinates	×
215	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
216	my = quadPointsSphere[idx_quad][0];
217	phi = quadPointsSphere[idx_quad][1];	×
218	moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
219	}	×
220	delete tempBase;	×
221		×
222	double s = 0.1;	×
223	double enterPositionX = 2.5; // 0.0;
224	double enterPositionY = 5.8;	×
225	double meanDir = M_PI / 2;
226	double s_ang = 0.1;	×
227		×
228	for( unsigned idx_cell = 0; idx_cell < cellMids.size(); ++idx_cell ) {	×
229	double x = cellMids[idx_cell][0] - enterPositionX;	×
230	double y = cellMids[idx_cell][1] - enterPositionY;	×
231	// anisotropic, forward-directed particle inflow
232	// for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {	×
233	// double theta = acos(sqrt(1-quadPointsSphere[idx_quad][0]quadPointsSphere[idx_quad][0])cos(quadPointsSphere[idx_quad][1]));	×
234	// if(theta > M_PI/3 && theta < 2*M_PI/3&&quadPointsSphere[idx_quad][0]<0) {	×
235	// cellKineticDensity[idx_quad] = std::max( 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ),
236	// epsilon );
237	// }
238	// }
239
240	// normal distribution also in angle
241
242	for( unsigned idx_quad = 0; idx_quad < _settings->GetNQuadPoints(); idx_quad++ ) {
243	if( quadPointsSphere[idx_quad][0] < 0 ) {
244	double theta =
245	acos( sqrt( 1 - quadPointsSphere[idx_quad][0] * quadPointsSphere[idx_quad][0] ) * cos( quadPointsSphere[idx_quad][1] ) );
246	double ang = theta - meanDir;	×
247	cellKineticDensity[idx_quad] =	×
248	std::max( 1.0 / ( s * s_ang * sqrt( 8 * M_PI * M_PI * M_PI ) ) * std::exp( -( x * x + y * y ) / ( 2 * s * s ) ) *
249	std::exp( -( ang * ang ) / ( 2 * s_ang ) ),	×
250	epsilon );	×
251	}	×
252	}	×
253		×
254	// Compute moments of this kinetic density	×
255	for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
256	initialSolution[idx_cell] += cellKineticDensity[idx_quad] * w[idx_quad] * moments[idx_quad];
257	}
258	}
259	delete quad;	×
260	return initialSolution;	×
261	}
262	}

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