snsolver.cpp 8.02 KB
Newer Older
steffen.schotthoefer's avatar
steffen.schotthoefer committed
1
#include "solvers/snsolver.h"
2
3
4
#include "common/config.h"
#include "common/io.h"
#include "common/mesh.h"
5
#include "fluxes/numericalflux.h"
6
#include "kernels/scatteringkernelbase.h"
7
#include "quadratures/quadraturebase.h"
8

9
10
// externals
#include "spdlog/spdlog.h"
11
#include <mpi.h>
Jonas Kusch's avatar
Jonas Kusch committed
12

13
14
SNSolver::SNSolver( Config* settings ) : Solver( settings ) {

15
16
17
18
19
    _quadPoints = _quadrature->GetPoints();
    _weights    = _quadrature->GetWeights();

    ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );
    _scatteringKernel   = k->GetScatteringKernel();
20
21
    delete k;
}
Jonas Kusch's avatar
Jonas Kusch committed
22

23
void SNSolver::Solve() {
24
    auto log = spdlog::get( "event" );
25
26

    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
27
    VectorVector psiNew = _sol;
vavrines's avatar
vavrines committed
28

vavrines's avatar
vavrines committed
29
30
31
32
33
34
35
    // reconstruction order
    unsigned reconsOrder = _settings->GetReconsOrder();

    // left and right angular flux of interface, used in numerical flux evaluation
    double psiL;
    double psiR;

vavrines's avatar
vavrines committed
36
    // derivatives of angular flux in x and y directions
vavrines's avatar
vavrines committed
37
38
    VectorVector psiDx( _nCells, Vector( _nq, 0.0 ) );
    VectorVector psiDy( _nCells, Vector( _nq, 0.0 ) );
vavrines's avatar
vavrines committed
39

vavrines's avatar
vavrines committed
40
    // geometric variables for derivatives computation
vavrines's avatar
vavrines committed
41
42
43
44
45
    auto nodes         = _mesh->GetNodes();
    auto cells         = _mesh->GetCells();
    auto cellMidPoints = _mesh->GetCellMidPoints();

    // center location of cell interfaces
vavrines's avatar
vavrines committed
46
    std::vector<std::vector<Vector>> interfaceMidPoints( _nCells, std::vector<Vector>( _mesh->GetNumNodesPerCell(), Vector( 2, 1e-10 ) ) );
vavrines's avatar
vavrines committed
47
48
    for( unsigned i = 0; i < _nCells; ++i ) {
        for( unsigned k = 0; k < _mesh->GetDim(); ++k ) {
vavrines's avatar
vavrines committed
49
50
            for( unsigned j = 0; j < _neighbors[i].size() - 1; ++j ) {
                interfaceMidPoints[i][j][k] = 0.5 * ( nodes[cells[i][j]][k] + nodes[cells[i][j + 1]][k] );
vavrines's avatar
vavrines committed
51
            }
vavrines's avatar
vavrines committed
52
            interfaceMidPoints[i][_neighbors[i].size() - 1][k] = 0.5 * ( nodes[cells[i][_neighbors[i].size() - 1]][k] + nodes[cells[i][0]][k] );
vavrines's avatar
vavrines committed
53
54
55
56
57
58
59
60
61
62
63
64
65
        }
    }

    // distance between cell center to interface center
    VectorVector cellDx( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
    VectorVector cellDy( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
    for( unsigned i = 0; i < _nCells; ++i ) {
        for( unsigned j = 0; j < _mesh->GetNumNodesPerCell(); ++j ) {
            cellDx[i][j] = interfaceMidPoints[i][j][0] - cellMidPoints[i][0];
            cellDy[i][j] = interfaceMidPoints[i][j][1] - cellMidPoints[j][1];
        }
    }

vavrines's avatar
vavrines committed
66
    double dFlux = 1e10;
67
68
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );
69
70
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
71
    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux" );
Jonas Kusch's avatar
Jonas Kusch committed
72

73
    // loop over energies (pseudo-time)
74
    for( unsigned n = 0; n < _nEnergies; ++n ) {
vavrines's avatar
vavrines committed
75
76
77

        // reconstruct slopes for higher order solver
        if( reconsOrder > 1 ) {
78
            _mesh->ReconstructSlopesU( _nq, psiDx, psiDy, _sol );    // unstructured reconstruction
vavrines's avatar
vavrines committed
79
80
            //_mesh->ReconstructSlopesS( _nq, psiDx, psiDy, _psi );    // structured reconstruction (not stable currently)
        }
vavrines's avatar
vavrines committed
81

82
        // loop over all spatial cells
Jonas Kusch's avatar
Jonas Kusch committed
83
        for( unsigned j = 0; j < _nCells; ++j ) {
84
            if( _boundaryCells[j] == BOUNDARY_TYPE::DIRICHLET ) continue;
85
86
            // loop over all ordinates
            for( unsigned k = 0; k < _nq; ++k ) {
87
                psiNew[j][k] = 0.0;
88
89
                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
                for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
90
                    // store flux contribution on psiNew_sigmaS to save memory
91
                    if( _boundaryCells[j] == BOUNDARY_TYPE::NEUMANN && _neighbors[j][l] == _nCells )
92
                        psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[j][k], _normals[j][l] );
vavrines's avatar
vavrines committed
93
                    else {
vavrines's avatar
vavrines committed
94
95
                        switch( reconsOrder ) {
                            // first order solver
96
                            case 1: psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[_neighbors[j][l]][k], _normals[j][l] ); break;
vavrines's avatar
vavrines committed
97
98
99
                            // second order solver
                            case 2:
                                // left status of interface
100
                                psiL = _sol[j][k] + psiDx[j][k] * ( interfaceMidPoints[j][l][0] - cellMidPoints[j][0] ) +
vavrines's avatar
vavrines committed
101
102
                                       psiDy[j][k] * ( interfaceMidPoints[j][l][1] - cellMidPoints[j][1] );
                                // right status of interface
103
                                psiR = _sol[_neighbors[j][l]][k] +
vavrines's avatar
vavrines committed
104
105
106
107
                                       psiDx[_neighbors[j][l]][k] * ( interfaceMidPoints[j][l][0] - cellMidPoints[_neighbors[j][l]][0] ) +
                                       psiDy[_neighbors[j][l]][k] * ( interfaceMidPoints[j][l][1] - cellMidPoints[_neighbors[j][l]][1] );
                                // positivity check (if not satisfied, deduce to first order)
                                if( psiL < 0.0 || psiR < 0.0 ) {
108
109
                                    psiL = _sol[j][k];
                                    psiR = _sol[_neighbors[j][l]][k];
vavrines's avatar
vavrines committed
110
111
112
113
114
115
116
                                }
                                // flux evaluation
                                psiNew[j][k] += _g->Flux( _quadPoints[k], psiL, psiR, _normals[j][l] );
                                break;
                            // higher order solver
                            case 3: std::cout << "higher order is WIP" << std::endl; break;
                            // default: first order solver
117
                            default: psiNew[j][k] += _g->Flux( _quadPoints[k], _sol[j][k], _sol[_neighbors[j][l]][k], _normals[j][l] );
vavrines's avatar
vavrines committed
118
119
                        }
                    }
Jonas Kusch's avatar
Jonas Kusch committed
120
121
                }
                // time update angular flux with numerical flux and total scattering cross section
122
                psiNew[j][k] = _sol[j][k] - ( _dE / _areas[j] ) * psiNew[j][k] - _dE * _sigmaT[n][j] * _sol[j][k];
Jonas Kusch's avatar
Jonas Kusch committed
123
124
            }
            // compute scattering effects
125
            psiNew[j] += _dE * _sigmaS[n][j] * _scatteringKernel * _sol[j];    // multiply scattering matrix with psi
126

127
            // TODO: figure out a more elegant way
128
            // add external source contribution
129
130
131
132
133
134
135
136
137
138
139
140
            if( _Q.size() == 1u ) {            // constant source for all energies
                if( _Q[0][j].size() == 1u )    // isotropic source
                    psiNew[j] += _dE * _Q[0][j][0];
                else
                    psiNew[j] += _dE * _Q[0][j];
            }
            else {
                if( _Q[0][j].size() == 1u )    // isotropic source
                    psiNew[j] += _dE * _Q[n][j][0];
                else
                    psiNew[j] += _dE * _Q[n][j];
            }
Jonas Kusch's avatar
Jonas Kusch committed
141
        }
142
        _sol = psiNew;
143
        for( unsigned i = 0; i < _nCells; ++i ) {
144
            fluxNew[i]       = dot( _sol[i], _weights );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
145
            _solverOutput[i] = fluxNew[i];
146
        }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
147
        // Save( n );
148
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
149
        fluxOld = fluxNew;
150
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[n], dFlux );
Jonas Kusch's avatar
Jonas Kusch committed
151
152
153
    }
}

154
void SNSolver::Save() const {
155
    std::vector<std::string> fieldNames{ "flux" };
156
157
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };

158
159
    std::vector<double> flux( _nCells, 0.0 );
    for( unsigned i = 0; i < _nCells; ++i ) {
160
        flux[i] = dot( _sol[i], _weights );
161
162
    }
    std::vector<std::vector<double>> scalarField( 1, flux );
163
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
164
    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
165
}
steffen.schotthoefer's avatar
steffen.schotthoefer committed
166
167

void SNSolver::Save( int currEnergy ) const {
168
    std::vector<std::string> fieldNames{ "flux" };
169
170
    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };

steffen.schotthoefer's avatar
steffen.schotthoefer committed
171
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
172
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
173
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNamesWrapper, _mesh );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
174
}