csdpnsolver.cpp 9.94 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"
Steffen Schotthöfer's avatar
update    
Steffen Schotthöfer committed
11
#include "toolboxes/textprocessingtoolbox.h"
jannick.wolters's avatar
jannick.wolters committed
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
    //
jannick.wolters's avatar
jannick.wolters committed
36
37
38
39
40
    double minE = 5e-5;
    double maxE = 1.0;
    _dE         = ComputeTimeStep( _settings->GetCFL() );
    _nEnergies  = std::ceil( ( maxE - minE ) / _dE );
    _energies   = blaze::linspace( _nEnergies, maxE, minE );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
41
42
43
44
45
    //_angle    = _energies;
    //
    //_sigmaSE  = { Matrix( _nEnergies, 0.0 ) };
    //_sigmaTE  = Vector( _nEnergies, 0.0 );
    //
jannick.wolters's avatar
jannick.wolters committed
46
47

    Vector pos_beam = Vector{ 0.5, 0.5 };
jannick.wolters's avatar
jannick.wolters committed
48
    VectorVector IC( _nCells, Vector( _nSystem ) );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
49
50
51
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        double x            = _cellMidPoints[idx_cell][0];
        double y            = _cellMidPoints[idx_cell][1];
52
        const double stddev = .005;
jannick.wolters's avatar
jannick.wolters committed
53
        double f            = normpdf( x, pos_beam[0], stddev ) * normpdf( y, pos_beam[1], stddev );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
54
55
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
            IC[idx_cell][idx_sys] = f * StarMAPmoments[idx_sys];    // must be VectorVector
56
        }
57
    }
58
59
    _sol    = IC;
    _solNew = IC;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
60

jannick.wolters's avatar
jannick.wolters committed
61
    _sigmaS = VectorVector( _nEnergies, Vector( _polyDegreeBasis ) );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
62
    for( unsigned idx_degree = 0; idx_degree < _polyDegreeBasis; ++idx_degree ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
63
        Vector xs_m = blaze::column( sigma_ref, idx_degree );    // Scattering cross section Moments
64
        Interpolation interp( E_ref, xs_m );
jannick.wolters's avatar
jannick.wolters committed
65
66
67
68
69
70
71
72
73
74
75
        auto tmp = interp( _energies );
        for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
            _sigmaS[idx_energy][idx_degree] = tmp[idx_energy];
        }
    }

    _sigmaT = VectorVector( _nEnergies, Vector( _polyDegreeBasis ) );
    for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
        for( unsigned idx_degree = 0; idx_degree < _polyDegreeBasis; ++idx_degree ) {
            _sigmaT[idx_energy][idx_degree] = _sigmaS[idx_energy][0] - _sigmaS[idx_energy][idx_degree];
        }
76
    }
jannick.wolters's avatar
jannick.wolters committed
77

jannick.wolters's avatar
jannick.wolters committed
78
    Interpolation interpS( E_tab, S_tab );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
79
    _s = interpS( _energies );
80
81
82
83
84

    // SanityChecks
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        _density[idx_cell] = 1.0;
    }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
85
}
86

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
87
88
89
90
void CSDPNSolver::SolverPreprocessing() {
    // TODO
}

91
void CSDPNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
92
    // TODO
jannick.wolters's avatar
jannick.wolters committed
93
94
}

95
void CSDPNSolver::IterPostprocessing( unsigned idx_iter ) {
96
97
98
99
100
101
    // --- Update Solution ---
    _sol = _solNew;

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

102
    unsigned n = idx_iter;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
103
104
105
106
107
108
109
110
111
112
113
114
    // -- 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
    }
}
jannick.wolters's avatar
jannick.wolters committed
115

116
void CSDPNSolver::FluxUpdate() {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
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
144
145
146
    // _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
147
148
                                               _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
149
150
151
152
                                               _normals[idx_cell][idx_neighbor] );
            }
        }
    }
153
154
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
155
156
157
void CSDPNSolver::FVMUpdate( unsigned idx_energy ) {
// loop over all spatial cells
#pragma omp parallel for
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
158
159
160
161
    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
162
163
164
165
166
167
168
169
        // for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                int idx_sys = GlobalIndex( idx_l, idx_k );

                _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] -
                                             ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
                                             - _dE * _sol[idx_cell][idx_sys] *
170
171
                                                   ( _sigmaT[idx_energy][idx_l] /* absorbtion influence */
                                                   );                           /* scattering influence */
172
            }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
173
        }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
174
        // Source Term
175
        // _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
176
177
    }
}
178

jannick.wolters's avatar
jannick.wolters committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
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
195

jannick.wolters's avatar
jannick.wolters committed
196
197
198
199
200
                _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
201
202
203
204
                _outputFields[idx_group].resize( 2 );
                _outputFieldNames[idx_group].resize( 2 );

                // Dose
jannick.wolters's avatar
jannick.wolters committed
205
                _outputFields[idx_group][0].resize( _nCells );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
206
207
208
209
                _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
210
211
                break;

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
212
            default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;
jannick.wolters's avatar
jannick.wolters committed
213
214
215
216
217
218
        }
    }
}

void CSDPNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
219
    double maxDose;
jannick.wolters's avatar
jannick.wolters committed
220
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
221
        ( idx_pseudoTime == _maxIter - 1 ) /* need sol at last iteration */ ) {
jannick.wolters's avatar
jannick.wolters committed
222
223
224
225
226

        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
227
                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
jannick.wolters's avatar
jannick.wolters committed
228
229
230
231
                    }
                    break;

                case MEDICAL:
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
232
233
234
235
236
237
238
239
240
                    // 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
241
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
242
                        _outputFields[idx_group][1][idx_cell] /= maxDose;
jannick.wolters's avatar
jannick.wolters committed
243
244
245
                    }
                    break;

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
246
                default: ErrorMessages::Error( "Volume Output Group not defined for CSD_SN_FP_TRAFO Solver!", CURRENT_FUNCTION ); break;
jannick.wolters's avatar
jannick.wolters committed
247
248
249
250
            }
        }
    }
}