csdsnsolverfp.cpp 8.93 KB
Newer Older
jonas.kusch's avatar
jonas.kusch committed
1
2
3
#include "solvers/csdsnsolverfp.h"
#include "common/config.h"
#include "common/io.h"
jannick.wolters's avatar
jannick.wolters committed
4
#include "common/mesh.h"
jonas.kusch's avatar
jonas.kusch committed
5
6
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
7
#include "problems/icru.h"
jonas.kusch's avatar
jonas.kusch committed
8
9
10
11
12
13
14
15
16
17
18
#include "problems/problembase.h"
#include "quadratures/quadraturebase.h"

// externals
#include "spdlog/spdlog.h"
#include <mpi.h>

CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );

    // Set angle and energies
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
19
20
    _angle     = Vector( _settings->GetNQuadPoints(), 0.0 );    // my
    _energies  = Vector( _nEnergies, 0.0 );                     // equidistant
21
22
23
    _energyMin = 1e-4;
    //_energyMax = 10.0;
    _energyMax = 5.0;
jonas.kusch's avatar
jonas.kusch committed
24
25
26
    // write equidistant energy grid

    _dE        = ComputeTimeStep( settings->GetCFL() );
jonas.kusch's avatar
jonas.kusch committed
27
    _nEnergies = unsigned( ( _energyMax - _energyMin ) / _dE );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
28
    // std::cout << "nEnergies = " << _nEnergies << std::endl;
jonas.kusch's avatar
jonas.kusch committed
29
30
    _energies.resize( _nEnergies );
    for( unsigned n = 0; n < _nEnergies; ++n ) {
jonas.kusch's avatar
jonas.kusch committed
31
        _energies[n] = _energyMin + ( _energyMax - _energyMin ) / ( _nEnergies - 1 ) * n;
jonas.kusch's avatar
jonas.kusch committed
32
33
34
    }

    // create 1D quadrature
jannick.wolters's avatar
jannick.wolters committed
35
    unsigned nq            = _settings->GetNQuadPoints();
36
    QuadratureBase* quad1D = QuadratureBase::Create( QUAD_GaussLegendre1D, nq );
jannick.wolters's avatar
jannick.wolters committed
37
38
39
40
    Vector w               = quad1D->GetWeights();
    VectorVector muVec     = quad1D->GetPoints();
    Vector mu( nq );
    for( unsigned k = 0; k < nq; ++k ) {
jonas.kusch's avatar
jonas.kusch committed
41
42
43
44
        mu[k] = muVec[k][0];
    }

    // setup Laplace Beltrami matrix L in slab geometry
jannick.wolters's avatar
jannick.wolters committed
45
    _L = Matrix( nq, nq, 0.0 );
jonas.kusch's avatar
jonas.kusch committed
46
47

    double DMinus = 0.0;
jannick.wolters's avatar
jannick.wolters committed
48
49
    double DPlus  = 0.0;
    for( unsigned k = 0; k < nq; ++k ) {
jonas.kusch's avatar
jonas.kusch committed
50
        DMinus = DPlus;
jannick.wolters's avatar
jannick.wolters committed
51
52
53
54
        DPlus  = DMinus - 2 * mu[k] * w[k];
        if( k > 0 ) {
            _L( k, k - 1 ) = DMinus / ( mu[k] - mu[k - 1] ) / w[k];
            _L( k, k )     = -DMinus / ( mu[k] - mu[k - 1] ) / w[k];
jonas.kusch's avatar
jonas.kusch committed
55
        }
jannick.wolters's avatar
jannick.wolters committed
56
57
58
        if( k < nq - 1 ) {
            _L( k, k + 1 ) = DPlus / ( mu[k + 1] - mu[k] ) / w[k];
            _L( k, k ) += -DPlus / ( mu[k + 1] - mu[k] ) / w[k];
jonas.kusch's avatar
jonas.kusch committed
59
60
61
62
        }
    }

    // determine momente of Heney-Greenstein
63
    _xi = Matrix( 3, _nEnergies );
jonas.kusch's avatar
jonas.kusch committed
64

jonas.kusch's avatar
jonas.kusch committed
65
    _s = Vector( _nEnergies, 1.0 );
jonas.kusch's avatar
jonas.kusch committed
66

jonas.kusch's avatar
jonas.kusch committed
67
    _RT = true;
jonas.kusch's avatar
jonas.kusch committed
68

jonas.kusch's avatar
jonas.kusch committed
69
    // read in medical data if radiation therapy option selected
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
70
71
    if( _RT ) {
        ICRU database( abs( mu ), _energies );
jonas.kusch's avatar
jonas.kusch committed
72
73
74
        database.GetTransportCoefficients( _xi );
        database.GetStoppingPower( _s );
    }
jonas.kusch's avatar
jonas.kusch committed
75
76
77
78
79
80
81
82
83

    // recompute scattering kernel. TODO: add this to kernel function
    for( unsigned p = 0; p < _nq; ++p ) {
        for( unsigned q = 0; q < _nq; ++q ) {
            _scatteringKernel( p, q ) = 0.0;
        }
        _scatteringKernel( p, p ) = _weights[p];
    }

84
    _density = std::vector<double>( _nCells, 1.0 );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
85
    // exit(EXIT_SUCCESS);
jonas.kusch's avatar
jonas.kusch committed
86
87
88
}

void CSDSNSolverFP::Solve() {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
89
    // std::cout << "Solve Fokker-Planck with Heney-Greenstein kernel using " << _nEnergies << " energies " << std::endl;
jonas.kusch's avatar
jonas.kusch committed
90

jannick.wolters's avatar
jannick.wolters committed
91
    auto log      = spdlog::get( "event" );
jonas.kusch's avatar
jonas.kusch committed
92
93
    auto cellMids = _mesh->GetCellMidPoints();

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
94
    _density = _problem->GetDensity( cellMids );
95

jonas.kusch's avatar
jonas.kusch committed
96
    // setup IC and incoming BC on left
jannick.wolters's avatar
jannick.wolters committed
97
    // auto cellMids = _settings->GetCellMidPoints();
jonas.kusch's avatar
jonas.kusch committed
98
    _sol = std::vector<Vector>( _nCells, Vector( _nq, 0.0 ) );
jannick.wolters's avatar
jannick.wolters committed
99
    for( unsigned k = 0; k < _nq; ++k ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
100
        if( _quadPoints[k][0] > 0 && !_RT ) _sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
jonas.kusch's avatar
jonas.kusch committed
101
    }
jannick.wolters's avatar
jannick.wolters committed
102
103
    _boundaryCells[0]           = BOUNDARY_TYPE::DIRICHLET;
    _boundaryCells[_nCells - 1] = BOUNDARY_TYPE::DIRICHLET;
jonas.kusch's avatar
jonas.kusch committed
104

jonas.kusch's avatar
jonas.kusch committed
105
106
107
    Matrix identity( _nq, _nq, 0.0 );
    for( unsigned k = 0; k < _nq; ++k ) identity( k, k ) = 1.0;

jonas.kusch's avatar
jonas.kusch committed
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
    VectorVector psiNew( _nCells, Vector( _nq, 0.0 ) );
    double dFlux = 1e10;
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );

    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    if( rank == 0 ) log->info( "{:10}   {:10}", "E", "dFlux" );

    // revert energies and stopping power
    Vector sSave = _s;
    for( unsigned n = 0; n < _nEnergies; ++n ) {
        _energies[n] = _energies[_nEnergies - 1] - _energies[n];
        _s[n]        = sSave[_nEnergies - 1 - n];
    }

    for( unsigned j = 0; j < _nCells; ++j ) {
        for( unsigned k = 0; k < _nq; ++k ) {
            fluxOld[j] += _weights[k] * _sol[j][k] * _s[0];
        }
    }

    // determine minimal density for CFL computation
    double densityMin = _density[0];
    for( unsigned j = 1; j < _nCells; ++j ) {
        if( densityMin > _density[j] ) densityMin = _density[j];
    }

    VectorVector psi1 = _sol;

    // cross sections do not need to be transformed to ETilde energy grid since e.g. TildeSigmaT(ETilde) = SigmaT(E(ETilde))

    // loop over energies (pseudo-time)
    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
144
        // set alpha and beta
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
145
146
        _alpha = _xi( 1, _nEnergies - n - 1 ) / 2.0 + _xi( 2, _nEnergies - n - 1 ) / 8.0;
        _beta  = _xi( 2, _nEnergies - n - 1 ) / 8.0 / _xi( 1, _nEnergies - n - 1 );
jonas.kusch's avatar
jonas.kusch committed
147
148
149
150

        _IL = identity - _beta * _L;

        // write BC
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
151
        if( _RT ) {
jonas.kusch's avatar
jonas.kusch committed
152
153
            for( unsigned k = 0; k < _nq; ++k ) {
                if( _quadPoints[k][0] > 0 ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
154
                    _sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp( -50 * pow( _energyMax - _energies[n], 2 ) );
jonas.kusch's avatar
jonas.kusch committed
155
156
157
158
                }
            }
        }

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
159
160
161
        // loop over all spatial cells
        // scattering step
        //#pragma omp parallel for
jonas.kusch's avatar
jonas.kusch committed
162
163
        for( unsigned j = 0; j < _nCells; ++j ) {
            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
jannick.wolters's avatar
jannick.wolters committed
164
165
            psi1[j] = blaze::solve( _IL, _sol[j] );
            psi1[j] = _alpha * _L * psi1[j];
jonas.kusch's avatar
jonas.kusch committed
166
        }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
167
168
        // advection step
        //#pragma omp parallel for
jonas.kusch's avatar
jonas.kusch committed
169
170
171
172
173
174
175
176
177
178
179
180
        for( unsigned j = 0; j < _nCells; ++j ) {
            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
            // loop over all ordinates
            for( unsigned i = 0; i < _nq; ++i ) {
                psiNew[j][i] = 0.0;
                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
                for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[j].size(); ++idx_neighbor ) {
                    // store flux contribution on psiNew_sigmaS to save memory
                    if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][idx_neighbor] == _nCells )
                        continue;    // adiabatic wall, add nothing
                    else
                        psiNew[j][i] += _g->Flux( _quadPoints[i], _sol[j][i], _sol[_neighbors[j][idx_neighbor]][i], _normals[j][idx_neighbor] ) /
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
181
                                        _areas[j] / ( _density[j] * _s[n + 1] );
jonas.kusch's avatar
jonas.kusch committed
182
183
                }
                // time update angular flux with numerical flux and total scattering cross section
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
184
185
                psiNew[j][i] =
                    _sol[j][i] - _dE * psiNew[j][i] + _dE * psi1[j][i] / ( _density[j] * _s[n + 1] ) + ( _s[n] / _s[n + 1] - 1 ) * _sol[j][i];
jonas.kusch's avatar
jonas.kusch committed
186
187
            }
        }
jonas.kusch's avatar
jonas.kusch committed
188
189
190
191
        for( unsigned j = 0; j < _nCells; ++j ) {
            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
            _sol[j] = psiNew[j];
        }
jonas.kusch's avatar
jonas.kusch committed
192

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
193
        //#pragma omp parallel for
jonas.kusch's avatar
jonas.kusch committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
        for( unsigned j = 0; j < _nCells; ++j ) {
            fluxNew[j] = 0.0;
            for( unsigned k = 0; k < _nq; ++k ) {
                fluxNew[j] += _weights[k] * _sol[j][k] * _s[n + 1];
            }
            _dose[j] += 0.5 * _dE * ( fluxNew[j] + fluxOld[j] );    // update dose with trapezoidal rule
            _solverOutput[j] = fluxNew[j];
        }

        // Save( n );
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
        fluxOld = fluxNew;
        if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}  {:01.5e}", _energies[n], _dE / densityMin, dFlux );
        if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
    }
    Save( 1 );
}

void CSDSNSolverFP::Save() const {
    std::vector<std::string> fieldNames{ "dose", "normalized dose" };
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
214
215
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };

jonas.kusch's avatar
jonas.kusch committed
216
217
218
219
220
    std::vector<std::vector<double>> dose( 1, _dose );
    std::vector<std::vector<double>> normalizedDose( 1, _dose );
    double maxDose = *std::max_element( _dose.begin(), _dose.end() );
    for( unsigned i = 0; i < _dose.size(); ++i ) normalizedDose[0][i] /= maxDose;
    std::vector<std::vector<std::vector<double>>> results{ dose, normalizedDose };
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
221
    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
jonas.kusch's avatar
jonas.kusch committed
222
223
224
225
}

void CSDSNSolverFP::Save( int currEnergy ) const {
    std::vector<std::string> fieldNames{ "flux" };
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
226
227
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };

jonas.kusch's avatar
jonas.kusch committed
228
229
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
230
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNamesWrapper, _mesh );
jonas.kusch's avatar
jonas.kusch committed
231
}