csdsnsolver.cpp 7.3 KB
Newer Older
1
#include "solvers/csdsnsolver.h"
2
3
#include "common/config.h"
#include "common/io.h"
4
#include "fluxes/numericalflux.h"
5
#include "kernels/scatteringkernelbase.h"
6
#include "problems/problembase.h"
7

8
9
10
11
// externals
#include "spdlog/spdlog.h"
#include <mpi.h>

12
13
14
15
CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );

    // Set angle and energies
Jonas Kusch's avatar
Jonas Kusch committed
16
17
    _angle           = Vector( _settings->GetNQuadPoints(), 0.0 );    // my
    _energies        = Vector( _nEnergies, 0.0 );                     // equidistant
18
    double energyMin = 1e-1;
Jonas Kusch's avatar
Jonas Kusch committed
19
    double energyMax = 5e0;
Jonas Kusch's avatar
Jonas Kusch committed
20
    // write equidistant energy grid
21
22
23
24

    _dE        = ComputeTimeStep( settings->GetCFL() );
    _nEnergies = unsigned( ( energyMax - energyMin ) / _dE );
    _energies.resize( _nEnergies );
Jonas Kusch's avatar
Jonas Kusch committed
25
26
27
28
29
30
31
    for( unsigned n = 0; n < _nEnergies; ++n ) {
        _energies[n] = energyMin + ( energyMax - energyMin ) / ( _nEnergies - 1 ) * n;
    }
    // write mu grid
    for( unsigned k = 0; k < _settings->GetNQuadPoints(); ++k ) {
        _angle[k] = _quadPoints[k][2];
    }
32

33
34
    _sigmaSE = _problem->GetScatteringXSE( _energies, _angle );
    _sigmaTE = _problem->GetTotalXSE( _energies );
Jonas Kusch's avatar
Jonas Kusch committed
35
    _s       = _problem->GetStoppingPower( _energies );
36
    _Q       = _problem->GetExternalSource( _energies );
Jonas Kusch's avatar
Jonas Kusch committed
37

38
    // Get patient density
jannick.wolters's avatar
jannick.wolters committed
39
    _density = Vector( _nCells, 1.0 );
40
}
Jonas Kusch's avatar
Jonas Kusch committed
41
42
43
44
45

void CSDSNSolver::Solve() {
    auto log = spdlog::get( "event" );

    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
Jonas Kusch's avatar
Jonas Kusch committed
46
47
    VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
    double dFlux = 1e10;
Jonas Kusch's avatar
Jonas Kusch committed
48
49
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );
50
    for( unsigned j = 0; j < _nCells; ++j ) {
51
        fluxOld[j] = dot( _sol[j], _weights );
52
    }
Jonas Kusch's avatar
Jonas Kusch committed
53
54
55
56
57
58
59
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if( rank == 0 ) log->info( "{:10}   {:10}", "E", "dFlux" );

    // do substitution from psi to psiTildeHat (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
    for( unsigned j = 0; j < _nCells; ++j ) {
        for( unsigned k = 0; k < _nq; ++k ) {
60
            _sol[j][k] = _sol[j][k] * _density[j] * _s[_nEnergies - 1];    // note that _s[_nEnergies - 1] is stopping power at highest energy
Jonas Kusch's avatar
Jonas Kusch committed
61
62
63
64
65
        }
    }

    // store transformed energies ETilde instead of E in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
    double tmp = 0.0;
Jonas Kusch's avatar
Jonas Kusch committed
66
67
68
    for( unsigned n = 0; n < _nEnergies; ++n ) {
        tmp          = tmp + _dE / _s[n];
        _energies[n] = tmp;
Jonas Kusch's avatar
Jonas Kusch committed
69
70
71
    }

    // store transformed energies ETildeTilde instead of ETilde in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
Jonas Kusch's avatar
Jonas Kusch committed
72
73
    for( unsigned n = 0; n < _nEnergies; ++n ) {
        _energies[n] = _energies[_nEnergies - 1] - _energies[n];
Jonas Kusch's avatar
Jonas Kusch committed
74
75
    }

Jonas Kusch's avatar
Jonas Kusch committed
76
77
78
79
80
81
    // determine minimal density for CFL computation
    double densityMin = _density[0];
    for( unsigned j = 1; j < _nCells; ++j ) {
        if( densityMin > _density[j] ) densityMin = _density[j];
    }

Jonas Kusch's avatar
Jonas Kusch committed
82
83
    // cross sections do not need to be transformed to ETilde energy grid since e.g. TildeSigmaT(ETilde) = SigmaT(E(ETilde))

Jonas Kusch's avatar
Jonas Kusch committed
84
    // loop over energies (pseudo-time)
Jonas Kusch's avatar
Jonas Kusch committed
85
86
    for( unsigned n = 0; n < _nEnergies - 1; ++n ) {
        _dE = fabs( _energies[n + 1] - _energies[n] );    // is the sign correct here?
Jonas Kusch's avatar
Jonas Kusch committed
87
        // loop over all spatial cells
Jonas Kusch's avatar
Jonas Kusch committed
88
89
        for( unsigned j = 0; j < _nCells; ++j ) {
            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
Jonas Kusch's avatar
Jonas Kusch committed
90
            // loop over all ordinates
Jonas Kusch's avatar
Jonas Kusch committed
91
92
            for( unsigned i = 0; i < _nq; ++i ) {
                psiNew[j][i] = 0.0;
Jonas Kusch's avatar
Jonas Kusch committed
93
                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
Jonas Kusch's avatar
Jonas Kusch committed
94
                for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
Jonas Kusch's avatar
Jonas Kusch committed
95
                    // store flux contribution on psiNew_sigmaS to save memory
Jonas Kusch's avatar
Jonas Kusch committed
96
                    if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
97
                        continue;    // adiabatic wall, add nothing
Jonas Kusch's avatar
Jonas Kusch committed
98
                    else
Jonas Kusch's avatar
Jonas Kusch committed
99
100
101
102
                        psiNew[j][i] += _g->Flux( _quadPoints[i],
                                                  _sol[j][i] / _density[j],
                                                  _sol[_neighbors[j][idx_neighbor]][i] / _density[_neighbors[j][idx_neighbor]],
                                                  _normals[j][idx_neighbor] );
Jonas Kusch's avatar
Jonas Kusch committed
103
104
                }
                // time update angular flux with numerical flux and total scattering cross section
105
                psiNew[j][i] = _sol[j][i] - ( _dE / _areas[j] ) * psiNew[j][i] - _dE * _sigmaTE[n] * _sol[j][i];
Jonas Kusch's avatar
Jonas Kusch committed
106
            }
Jonas Kusch's avatar
Jonas Kusch committed
107
            // compute scattering effects (_scatteringKernel is simply multiplication with quad weights)
108
            psiNew[j] += _dE * ( blaze::trans( _sigmaSE[n] ) * _scatteringKernel * _sol[j] );    // multiply scattering matrix with psi
109
            // TODO: Check if _sigmaS^T*psi is correct
Jonas Kusch's avatar
Jonas Kusch committed
110
111
112

            // TODO: figure out a more elegant way
            // add external source contribution
113
            if( _Q.size() == 1u ) {            // constant source for all energies
Jonas Kusch's avatar
Jonas Kusch committed
114
115
                if( _Q[0][j].size() == 1u )    // isotropic source
                    psiNew[j] += _dE * _Q[0][j][0] * _s[_nEnergies - n - 1];
Jonas Kusch's avatar
Jonas Kusch committed
116
                else
Jonas Kusch's avatar
Jonas Kusch committed
117
                    psiNew[j] += _dE * _Q[0][j] * _s[_nEnergies - n - 1];
Jonas Kusch's avatar
Jonas Kusch committed
118
119
            }
            else {
Jonas Kusch's avatar
Jonas Kusch committed
120
121
                if( _Q[0][j].size() == 1u )    // isotropic source
                    psiNew[j] += _dE * _Q[n][j][0] * _s[_nEnergies - n - 1];
Jonas Kusch's avatar
Jonas Kusch committed
122
                else
Jonas Kusch's avatar
Jonas Kusch committed
123
                    psiNew[j] += _dE * _Q[n][j] * _s[_nEnergies - n - 1];
124
            }
Jonas Kusch's avatar
Jonas Kusch committed
125
        }
126
        _sol = psiNew;
Jonas Kusch's avatar
Jonas Kusch committed
127
        // do backsubstitution from psiTildeHat to psi (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
Jonas Kusch's avatar
Jonas Kusch committed
128
129
130
        for( unsigned j = 0; j < _nCells; ++j ) {
            for( unsigned i = 0; i < _nq; ++i ) {
                psiNew[j][i] = _sol[j][i] * _density[j] * _s[_nEnergies - n - 1];    // note that _s[0] is stopping power at lowest energy
Jonas Kusch's avatar
Jonas Kusch committed
131
132
133
            }
        }

134
135
136
137
        for( unsigned j = 0; j < _nCells; ++j ) {
            fluxNew[j] = dot( psiNew[j], _weights );
            _dose[j] += 0.5 * _dE * ( fluxNew[j] + fluxOld[j] ) / _density[j];    // update dose with trapezoidal rule
            _solverOutput[j] = fluxNew[j];
Jonas Kusch's avatar
Jonas Kusch committed
138
        }
139

jannick.wolters's avatar
jannick.wolters committed
140
        // Save( n );
Jonas Kusch's avatar
Jonas Kusch committed
141
142
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
        fluxOld = fluxNew;
Jonas Kusch's avatar
Jonas Kusch committed
143
        if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}  {:01.5e}", _energies[n], _dE / densityMin, dFlux );
144
        if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
Jonas Kusch's avatar
Jonas Kusch committed
145
146
147
148
    }
}

void CSDSNSolver::Save() const {
149
    std::vector<std::string> fieldNames{ "dose" };
150
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
151
    std::vector<std::vector<double>> scalarField( 1, _dose );
Jonas Kusch's avatar
Jonas Kusch committed
152
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
153
    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
Jonas Kusch's avatar
Jonas Kusch committed
154
155
156
}

void CSDSNSolver::Save( int currEnergy ) const {
157
158
    std::vector<std::string> fieldNames{ "dose" };
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
Jonas Kusch's avatar
Jonas Kusch committed
159
160
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
161
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNamesWrapper, _mesh );
Jonas Kusch's avatar
Jonas Kusch committed
162
}