csdpnsolver.cpp 8.78 KB
Newer Older
1
#include "solvers/csdpnsolver.h"
jannick.wolters's avatar
jannick.wolters committed
2
#include "common/config.h"
jannick.wolters's avatar
jannick.wolters committed
3
#include "common/globalconstants.h"
jannick.wolters's avatar
jannick.wolters committed
4
#include "common/io.h"
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
5
#include "common/mesh.h"
jannick.wolters's avatar
jannick.wolters committed
6
7
8
9
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
10
#include "solvers/csdpn_starmap_constants.h"
jannick.wolters's avatar
jannick.wolters committed
11
12
13
14
15

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

jannick.wolters's avatar
jannick.wolters committed
16
17
18
19
double normpdf( double x, double mu, double sigma ) {
    return INV_SQRT_2PI / sigma * std::exp( -( ( x - mu ) * ( x - mu ) ) / ( 2.0 * sigma * sigma ) );
}

20
21
22
23
24
25
26
27
28
29
30
Vector Time2Energy( const Vector& t, const double E_CutOff ) {
    Interpolation interp( E_trans, E_tab );
    Interpolation interp2( E_tab, E_trans );
    return blaze::max( 0, interp( interp2( E_CutOff, 0 ) - t ) );
}

Vector Energy2Time( const Vector& E, const double E_CutOff ) {
    Interpolation interp( E_tab, E_trans );
    return blaze::max( 0, interp( E_CutOff - E ) );
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
31
32
33
34
CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );

    // only dummies for compilation
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
35
36
37
38
39
40
41
    //
    //_energies = Vector( _nEnergies, 0.0 );
    //_angle    = _energies;
    //
    //_sigmaSE  = { Matrix( _nEnergies, 0.0 ) };
    //_sigmaTE  = Vector( _nEnergies, 0.0 );
    //
jannick.wolters's avatar
jannick.wolters committed
42
43

    Vector pos_beam = Vector{ 0.5, 0.5 };
jannick.wolters's avatar
jannick.wolters committed
44
    VectorVector IC( _nCells, Vector( _nSystem ) );
45
46
47
48
    for( unsigned i = 0; i < _nCells; ++i ) {
        double x            = _cellMidPoints[i][0];
        double y            = _cellMidPoints[i][1];
        const double stddev = .005;
jannick.wolters's avatar
jannick.wolters committed
49
        double f            = normpdf( x, pos_beam[0], stddev ) * normpdf( y, pos_beam[1], stddev );
50
        for( unsigned j = 0; j < _nSystem; j++ ) {
jannick.wolters's avatar
jannick.wolters committed
51
            IC[i][j] = f * StarMAPmoments[j];    // must be VectorVector
52
        }
53
    }
jannick.wolters's avatar
jannick.wolters committed
54

jannick.wolters's avatar
jannick.wolters committed
55
    Matrix sigma_t( _energies.size(), sigma_ref.rows() );
56
57
58
59
60
    for( unsigned i = 0; i < _nSystem; ++i ) {
        Vector xs_m = blaze::column( sigma_ref, i );
        Interpolation interp( E_ref, xs_m );
        // blaze::column( sigma_t, i ) = interp( _energies );
    }
jannick.wolters's avatar
jannick.wolters committed
61

jannick.wolters's avatar
jannick.wolters committed
62
    Interpolation interpS( E_tab, S_tab );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
63
    _s = interpS( _energies );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
64
}
65

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
66
67
68
69
void CSDPNSolver::SolverPreprocessing() {
    // TODO
}

70
void CSDPNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
71
    // TODO
jannick.wolters's avatar
jannick.wolters committed
72
73
}

74
void CSDPNSolver::IterPostprocessing( unsigned idx_iter ) {
75
76
77
78
79
80
    // --- Update Solution ---
    _sol = _solNew;

    // --- Compute Flux for solution and Screen Output ---
    ComputeRadFlux();

81
    unsigned n = idx_iter;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
82
83
84
85
86
87
88
89
90
91
92
93
94
    // -- Compute Dose
    for( unsigned j = 0; j < _nCells; ++j ) {
        if( n > 0 ) {
            _dose[j] += 0.5 * _dE * ( _fluxNew[j] * _s[_nEnergies - n - 1] + _flux[j] * _s[_nEnergies - n] ) /
                        _density[j];    // update dose with trapezoidal rule
        }
        else {
            _dose[j] += _dE * _fluxNew[j] * _s[_nEnergies - n - 1] / _density[j];
        }
        _solverOutput[j] = _fluxNew[j];    // Carefull here
        _flux[j]         = _fluxNew[j];    // Carefull here
    }
}
jannick.wolters's avatar
jannick.wolters committed
95

96
void CSDPNSolver::FluxUpdate() {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    // _mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol );

#pragma omp parallel for
    // Loop over all spatial cells
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {

        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;

        // Reset temporary variable psiNew
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
            _solNew[idx_cell][idx_sys] = 0.0;
        }

        // Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
        for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cell].size(); idx_neighbor++ ) {

            // Compute flux contribution and store in psiNew to save memory
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
                _solNew[idx_cell] += _g->Flux(
                    _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
            else {

                // first order solver
                _solNew[idx_cell] += _g->Flux( _AxPlus,
                                               _AxMinus,
                                               _AyPlus,
                                               _AyMinus,
                                               _AzPlus,
                                               _AzMinus,
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
127
128
                                               _sol[idx_cell] / _density[idx_cell],
                                               _sol[_neighbors[idx_cell][idx_neighbor]] / _density[_neighbors[idx_cell][idx_neighbor]],
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
129
130
131
132
                                               _normals[idx_cell][idx_neighbor] );
            }
        }
    }
133
134
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
135
136
137
void CSDPNSolver::FVMUpdate( unsigned idx_energy ) {
// loop over all spatial cells
#pragma omp parallel for
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
138
139
140
141
142
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // Flux update
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
143
            _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
144
145
146
                                         - _dE * _sol[idx_cell][idx_sys] *
                                               ( _sigmaT[idx_energy][idx_cell]                                 /* absorbtion influence */
                                                 + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
147
        }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
148
149
        // Source Term
        _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
150
151
    }
}
152

jannick.wolters's avatar
jannick.wolters committed
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
void CSDPNSolver::PrepareVolumeOutput() {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();

    _outputFieldNames.resize( nGroups );
    _outputFields.resize( nGroups );

    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
        // Prepare all Output Fields per group

        // Different procedure, depending on the Group...
        switch( _settings->GetVolumeOutput()[idx_group] ) {
            case MINIMAL:
                // Currently only one entry ==> rad flux
                _outputFields[idx_group].resize( 1 );
                _outputFieldNames[idx_group].resize( 1 );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
169

jannick.wolters's avatar
jannick.wolters committed
170
171
172
173
174
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
                break;

            case MEDICAL:
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
175
176
177
178
                _outputFields[idx_group].resize( 2 );
                _outputFieldNames[idx_group].resize( 2 );

                // Dose
jannick.wolters's avatar
jannick.wolters committed
179
                _outputFields[idx_group][0].resize( _nCells );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
180
181
182
183
                _outputFieldNames[idx_group][0] = "dose";
                // Normalized Dose
                _outputFields[idx_group][1].resize( _nCells );
                _outputFieldNames[idx_group][1] = "normalized dose";
jannick.wolters's avatar
jannick.wolters committed
184
185
                break;

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
186
            default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;
jannick.wolters's avatar
jannick.wolters committed
187
188
189
190
191
192
        }
    }
}

void CSDPNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
193
    double maxDose;
jannick.wolters's avatar
jannick.wolters committed
194
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
195
        ( idx_pseudoTime == _maxIter - 1 ) /* need sol at last iteration */ ) {
jannick.wolters's avatar
jannick.wolters committed
196
197
198
199
200

        for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
            switch( _settings->GetVolumeOutput()[idx_group] ) {
                case MINIMAL:
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
201
                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
jannick.wolters's avatar
jannick.wolters committed
202
203
204
205
                    }
                    break;

                case MEDICAL:
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
206
207
208
209
210
211
212
213
214
                    // Compute Dose
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        _outputFields[idx_group][0][idx_cell] += _dose[idx_cell];
                    }
                    // Compute normalized dose
                    _outputFields[idx_group][1] = _outputFields[idx_group][0];

                    maxDose = *std::max_element( _outputFields[idx_group][0].begin(), _outputFields[idx_group][0].end() );

jannick.wolters's avatar
jannick.wolters committed
215
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
216
                        _outputFields[idx_group][1][idx_cell] /= maxDose;
jannick.wolters's avatar
jannick.wolters committed
217
218
219
                    }
                    break;

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
220
                default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;
jannick.wolters's avatar
jannick.wolters committed
221
222
223
224
            }
        }
    }
}