csdpnsolver.cpp 8.89 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
CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
32
33
    blaze::transpose( sigma_ref );

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
34
35
36
    _dose = std::vector<double>( _settings->GetNCells(), 0.0 );

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

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

59
60
61
62
63
64
    Matrix sigma_t( _energies.size(), sigma_t.rows() );
    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
65

66
    Interpolation interpS( E_ref, S_tab );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
67
    _s = interpS( _energies );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
68
}
69

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
70
71
72
73
void CSDPNSolver::SolverPreprocessing() {
    // TODO
}

74
void CSDPNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
75
    // TODO
jannick.wolters's avatar
jannick.wolters committed
76
77
}

78
void CSDPNSolver::IterPostprocessing( unsigned idx_iter ) {
79
80
81
82
83
84
    // --- Update Solution ---
    _sol = _solNew;

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

85
    unsigned n = idx_iter;
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
86
87
88
89
90
91
92
93
94
95
96
97
98
    // -- 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
99

100
void CSDPNSolver::FluxUpdate() {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
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
127
128
129
130
    // _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
131
132
                                               _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
133
134
135
136
                                               _normals[idx_cell][idx_neighbor] );
            }
        }
    }
137
138
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
139
140
141
void CSDPNSolver::FVMUpdate( unsigned idx_energy ) {
// loop over all spatial cells
#pragma omp parallel for
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
142
143
144
145
146
    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
147
            _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
148
149
150
                                         - _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
151
        }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
152
153
        // Source Term
        _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
154
155
    }
}
156

jannick.wolters's avatar
jannick.wolters committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
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
173

jannick.wolters's avatar
jannick.wolters committed
174
175
176
177
178
                _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
179
180
181
182
                _outputFields[idx_group].resize( 2 );
                _outputFieldNames[idx_group].resize( 2 );

                // Dose
jannick.wolters's avatar
jannick.wolters committed
183
                _outputFields[idx_group][0].resize( _nCells );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
184
185
186
187
                _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
188
189
                break;

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

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

        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
205
                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
jannick.wolters's avatar
jannick.wolters committed
206
207
208
209
                    }
                    break;

                case MEDICAL:
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
210
211
212
213
214
215
216
217
218
                    // 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
219
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
220
                        _outputFields[idx_group][1][idx_cell] /= maxDose;
jannick.wolters's avatar
jannick.wolters committed
221
222
223
                    }
                    break;

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