csdsnsolver.cpp 6.47 KB
Newer Older
1
2
3
4
#include "fluxes/numericalflux.h"
#include "io.h"
#include "settings/config.h"

Jonas Kusch's avatar
Jonas Kusch committed
5
6
#include "solvers/csdsnsolver.h"

7
CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) { _dose = std::vector<double>( _settings->GetNCells(), 0.0 ); }
Jonas Kusch's avatar
Jonas Kusch committed
8
9
10
11
12
13
14
15
16

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?)
    VectorVector psiNew = _psi;
    double dFlux        = 1e10;
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );
17
18
19
    for( unsigned j = 0; j < _nCells; ++j ) {
        fluxOld[j] = dot( _psi[j], _weights );
    }
Jonas Kusch's avatar
Jonas Kusch committed
20
21
22
23
24
25
26
27
28
29
30
31
32
    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 ) {
            _psi[j][k] = _psi[j][k] * _density[j] * _s[_nEnergies - 1];    // note that _s[_nEnergies - 1] is stopping power at highest energy
        }
    }

    // store transformed energies ETilde instead of E in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
    double tmp = 0.0;
33
34
35
    for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
        tmp                   = tmp + _dE / _s[idx_energy];
        _energies[idx_energy] = tmp;
Jonas Kusch's avatar
Jonas Kusch committed
36
37
38
    }

    // store transformed energies ETildeTilde instead of ETilde in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
39
40
    for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
        _energies[idx_energy] = _energies[_nEnergies - 1] - _energies[idx_energy];
Jonas Kusch's avatar
Jonas Kusch committed
41
42
43
    }

    // loop over energies (pseudo-time)
44
45
    for( unsigned idx_energy = 1; idx_energy < _nEnergies; ++idx_energy ) {
        _dE = fabs( _energies[idx_energy] - _energies[idx_energy - 1] );    // is the sign correct here?
Jonas Kusch's avatar
Jonas Kusch committed
46
        // loop over all spatial cells
47
48
        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
Jonas Kusch's avatar
Jonas Kusch committed
49
            // loop over all ordinates
50
51
            for( unsigned idx_ord = 0; idx_ord < _nq; ++idx_ord ) {
                psiNew[idx_cell][idx_ord] = 0.0;
Jonas Kusch's avatar
Jonas Kusch committed
52
                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
53
                for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cell].size(); ++idx_neighbor ) {
Jonas Kusch's avatar
Jonas Kusch committed
54
                    // store flux contribution on psiNew_sigmaS to save memory
55
56
57
58
59
60
                    if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
                        continue;    // adiabatic wall, add nothing
                    // psiNew[idx_cell][idx_ord] += _g->Flux( _quadPoints[idx_ord] / _density[idx_cell],
                    //                                _psi[idx_cell][idx_ord],
                    //                              _psi[idx_cell][idx_ord],
                    //                            _normals[idx_cell][idx_neighbor] );
Jonas Kusch's avatar
Jonas Kusch committed
61
                    else
62
63
64
65
                        psiNew[idx_cell][idx_ord] += _g->Flux( _quadPoints[idx_ord] / _density[idx_cell],
                                                               _psi[idx_cell][idx_ord],
                                                               _psi[_neighbors[idx_cell][idx_neighbor]][idx_ord],
                                                               _normals[idx_cell][idx_neighbor] );
Jonas Kusch's avatar
Jonas Kusch committed
66
67
                }
                // time update angular flux with numerical flux and total scattering cross section
68
69
                psiNew[idx_cell][idx_ord] = _psi[idx_cell][idx_ord] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_ord] -
                                            _dE * _sigmaT[idx_energy][idx_cell] * _psi[idx_cell][idx_ord];
Jonas Kusch's avatar
Jonas Kusch committed
70
71
            }
            // compute scattering effects
72
            psiNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _psi[idx_cell];    // multiply scattering matrix with psi
Jonas Kusch's avatar
Jonas Kusch committed
73
74
75

            // TODO: figure out a more elegant way
            // add external source contribution
76
77
78
            if( _Q.size() == 1u ) {                   // constant source for all energies
                if( _Q[0][idx_cell].size() == 1u )    // isotropic source
                    psiNew[idx_cell] += _dE * _Q[0][idx_cell][0] * _s[_nEnergies - idx_energy - 1];
Jonas Kusch's avatar
Jonas Kusch committed
79
                else
80
                    psiNew[idx_cell] += _dE * _Q[0][idx_cell] * _s[_nEnergies - idx_energy - 1];
Jonas Kusch's avatar
Jonas Kusch committed
81
82
            }
            else {
83
84
                if( _Q[0][idx_cell].size() == 1u )    // isotropic source
                    psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0] * _s[_nEnergies - idx_energy - 1];
Jonas Kusch's avatar
Jonas Kusch committed
85
                else
86
                    psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell] * _s[_nEnergies - idx_energy - 1];
Jonas Kusch's avatar
Jonas Kusch committed
87
88
89
90
91
            }
        }
        _psi = psiNew;

        // do backsubstitution from psiTildeHat to psi (cf. Dissertation Kerstion Kuepper, Eq. 1.23)
92
93
94
95
        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
            for( unsigned idx_ord = 0; idx_ord < _nq; ++idx_ord ) {
                psiNew[idx_cell][idx_ord] = _psi[idx_cell][idx_ord] * _density[idx_cell] *
                                            _s[_nEnergies - idx_energy - 1];    // note that _s[0] is stopping power at lowest energy
Jonas Kusch's avatar
Jonas Kusch committed
96
97
98
            }
        }

99
100
101
102
        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
103
        }
104
105

        Save( idx_energy );
Jonas Kusch's avatar
Jonas Kusch committed
106
107
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
        fluxOld = fluxNew;
108
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], dFlux );
Jonas Kusch's avatar
Jonas Kusch committed
109
110
111
112
    }
}

void CSDSNSolver::Save() const {
113
114
    std::vector<std::string> fieldNames{ "dose" };
    std::vector<std::vector<double>> scalarField( 1, _dose );
Jonas Kusch's avatar
Jonas Kusch committed
115
116
117
118
119
120
121
122
123
124
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}

void CSDSNSolver::Save( int currEnergy ) const {
    std::vector<std::string> fieldNames{ "flux" };
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
}