Homec4science
Diffusion ConSol (master) Browse

ConSol/ConSol/ConSol/Reconstr/WENO_structured.cpp
759a8a6a6555master

Last Change
Raw File

WENO_structured.cpp
View Options

//
// Created by Fabian Moenkeberg on 2019-09-18.
//
#include "WENO_structured.hpp"
#include <assert.h>
#include "../Utility/MatrixVectorProd.hpp"
//#include "omp.h"
#include <omp.h>
#include <numeric>
#include "../Utility/ddGeneralizedNewton.hpp"
#include "../Utility/signPropFulfilled.hpp"
#include "../Utility/isWellDef.hpp"
#include "Utility/gauss-legendre-quadrature/gauss_legendre.hpp"
#include "H5Cpp.h"
#include "Utility/Hesthaven/stdWenoWeights.hpp"
void WENO_structured::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Configuration config){
this->model = model;
this->limiter = limiter;
this->interpBase = interpBase;
this->order = ceil((config.getOrder()-1)/2.0)+1;
this->nThreads = config.getNthreads();
int n = 1;
std::vector<double> w0(n/2 + (n%2));
std::vector<double> abs(n/2 + (n%2));
Qlambdas.resize(n, std::vector<double>(2));
QWeights.resize(n);
double *w = new double[n];
double *x = new double[n];
legendre_set(n,x,w);
if (order < 2){
Qlambdas[0][0] = 0.5; Qlambdas[0][1]= 0.5;
QWeights[0] = 1;
}else {
for (int i = 0; i < n; i++) {
Qlambdas[i][0] = (1 + x[i]) / 2;
Qlambdas[i][1] = (1 - x[i]) / 2;
QWeights[i] = w[i] / 2;
}
}
if (nThreads == 0){
nThreads = omp_get_max_threads();
}
std::cout << " Reconstruction: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n";
delete [] w;
delete [] x;
xCoord = mesh->getXCoord();
int ngc = mesh->getNgc();
nX = xCoord[0].size()-1-2*ngc; // Number of cells in x-dir.
nY = xCoord[1].size()-1-2*ngc; // Number of cells in y-dir.
// Initialize linear weights WENO reconstruction to get higher order by combining the low
// order reconstructions v_{i+1/2} = sum_{r=0}^{k-1} dw(r)v^{r}_{i+1/2}.
dw = LinearWeights(order,0);
Crec.resize(order + 1, std::vector<double> (order + 1));
for (int r = -1; r < order; r++){
Crec[r+1] = ReconstructWeights(order, r);
}
std::vector<std::vector<double>> beta0;
beta.resize(order, std::vector<std::vector<double>>(order, std::vector<double>(order,0)));
std::vector<double> xL;
for (int r = 0; r < order; r++){
xL.resize(order + 1);
for (int ii = 0; ii <= order; ii++){
xL[ii] = -0.5 + ii - r;
}
beta0 = betarcalc(xL, order);
for (int ii = 0; ii < order; ii ++){
for (int jj = 0; jj < order; jj ++){
beta[ii][jj][r] = beta0[ii][jj];
}
}
}
}
void WENO_structured::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U, MeshBase *mesh, double t){
int dimSyst = mesh->getNelem();
int nQuads = QWeights.size();
// std::vector<std::vector<std::vector<double>>> UReconstr0(nX,std::vector<std::vector<double>>(nY,std::vector<double>(dimSyst)));
omp_set_num_threads(nThreads);
int nNodes = 4;
std::vector<double> yloc;
// std::vector<double> rec(2);
std::vector<double> xloc;
// std::vector<double> uloc(2*(order-1) + 1);
#pragma omp parallel for
for (int i = 0; i < nX; i++) {
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int j = 0; j < order; j++){
// uloc.resize(order + j);
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(l+j)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(j+1)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i+(j)*nX][0][k] = rec[0];
UReconstr[0][i+(j+1)*nX][0][k] = rec[1];
}
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(nY+ngc - order - j + l)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(nY - j - 1)*nX][0][k] = U[i+ngc+((nY - j - 1)+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(nY - j)*nX][0][k] = U[i+ngc+((nY - j - 1)+ngc)*(nX+2*ngc)][k];
UReconstr[1][i+(nY - j - 1)*nX][0][k] = rec[0];
UReconstr[0][i+(nY - j)*nX][0][k] = rec[1];
}
}
for (int j = order; j < nY+1-order; j++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(l+j+ngc-1)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(j+1)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i+(j)*nX][0][k] = rec[0];
UReconstr[0][i+(j+1)*nX][0][k] = rec[1];
}
}
}
#pragma omp parallel for
for (int j = 0; j < nY; j++) {
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int i = 0; i < order; i++){
// uloc.resize(order + i);
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[(l+i)+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = rec[0];
UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = rec[1];
}
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[nX+ngc-order-i+l+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][nX-1 - i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][nX-1 - i+1 +(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][nX-1 - i+(j)*(nX+1) + nX*(nY+1)][0][k] = rec[0];
UReconstr[0][nX-1 - i+1 +(j)*(nX+1) + nX*(nY+1)][0][k] = rec[1];
}
}
for (int i = order; i < nX+1-order; i++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+l+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = rec[0];
UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = rec[1];
}
// if (i+(j)*nX + nX*(nY+1)> 119){
// std::cout << i+(j)*nX + nX*(nY+1) << "\n";
// }
}
}
mesh->updateEdges(UReconstr, t);
}
void WENO_structured::reconstruct(std::vector<std::vector<double>> &numFlux, NetFlux* numFluxFct, std::vector<std::vector<double>> &U, MeshBase *mesh, double t, double dt){
int dimSyst = mesh->getNelem();
std::vector<std::vector<std::vector<std::vector<double>>>> UReconstr(2,std::vector<std::vector<std::vector<double>>>(nY+1,std::vector<std::vector<double>>(1,std::vector<double>(dimSyst))));
omp_set_num_threads(nThreads);
std::vector<double> yloc;
std::vector<double> xloc;
double dX = xCoord[0][1]-xCoord[0][0];
double dY = xCoord[1][1]-xCoord[1][0];
std::vector<double> dxMin = mesh->getDxMin();
#pragma omp parallel for
for (int i = 0; i < nX; i++) {
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int j = 0; j < order; j++){
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(l+j)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(j+1)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][j][0][k] = rec[0];
UReconstr[0][j+1][0][k] = rec[1];
}
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(nY+ngc - order - j + l)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(nY - j - 1)*nX][0][k] = U[i+ngc+((nY - j - 1)+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(nY - j)*nX][0][k] = U[i+ngc+((nY - j - 1)+ngc)*(nX+2*ngc)][k];
UReconstr[1][(nY - j - 1)][0][k] = rec[0];
UReconstr[0][(nY - j)][0][k] = rec[1];
}
}
for (int j = order; j < nY+1-order; j++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+ngc+(l+j+ngc-1)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+(j+1)*nX][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][j][0][k] = rec[0];
UReconstr[0][j+1][0][k] = rec[1];
}
}
mesh->updateEdges(UReconstr, t, i);
mesh->updateEdges(UReconstr, t, i + nX);
for (int j = 0; j < nY; j++){
for (int k = 0; k < dimSyst; k++) {
numFlux[i+ngc+(j+ngc)*(nX+2*ngc)][k] = 0;
}
std::vector<double> flux1 = numFluxFct->Fevaly(UReconstr[0][j+1], UReconstr[1][j+1],t,dt, dxMin);
std::vector<double> flux2 = numFluxFct->Fevaly(UReconstr[0][j], UReconstr[1][j],t,dt, dxMin);
for (int k = 0; k < dimSyst; k++) {
numFlux[i+ngc+(j+ngc)*(nX+2*ngc)][k] -= (flux1[k]-flux2[k])/dY;
}
}
}
UReconstr.resize(2,std::vector<std::vector<std::vector<double>>>(nX+1,std::vector<std::vector<double>>(1,std::vector<double>(dimSyst))));
#pragma omp parallel for
for (int j = 0; j < nY; j++) {
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int i = 0; i < order; i++){
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[(l+i)+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = U[i+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i][0][k] = rec[0];
UReconstr[0][i+1][0][k] = rec[1];
}
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[nX+ngc-order-i+l+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][nX-1 - i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][nX-1 - i+1 +(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][nX-1 - i][0][k] = rec[0];
UReconstr[0][nX-1 - i+1][0][k] = rec[1];
}
}
for (int i = order; i < nX+1-order; i++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 0; l < 2*(order-1)+1; l++){
uloc[l] = U[i+l+(j+ngc)*(nX+2*ngc)][k];
}
WENO_ev(uloc, rec);
// UReconstr[1][i+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
// UReconstr[0][i+1+(j)*(nX+1) + nX*(nY+1)][0][k] = U[nX-i-1+ngc+(j+ngc)*(nX+2*ngc)][k];
UReconstr[1][i][0][k] = rec[0];
UReconstr[0][i+1][0][k] = rec[1];
}
}
mesh->updateEdges(UReconstr, t, j+2*nX);
mesh->updateEdges(UReconstr, t, j+2*nX + nY);
for (int i = 0; i < nX; i++){
std::vector<double> flux1 = numFluxFct->Fevalx(UReconstr[0][i+1], UReconstr[1][i+1],t,dt, dxMin);
std::vector<double> flux2 = numFluxFct->Fevalx(UReconstr[0][i], UReconstr[1][i],t,dt, dxMin);
for (int k = 0; k < dimSyst; k++) {
numFlux[i+ngc+(j+ngc)*(nX+2*ngc)][k] -= (flux1[k]-flux2[k])/dX;
}
}
}
}
//std::vector<int> WENO_structured::getStencil(int cellIndex, MeshBase *mesh, std::vector<std::vector<double>> *U, int k);{}
void WENO_structured::WENO_ev(std::vector<double> uloc, std::vector<double> &recloc){
double vareps = 1e-6;
// std::vector<double> upl(order,0);
// std::vector<double> uml(order,0);
// std::vector<double> betar00(order,0);
double upl;
double uml;
double betar;
double tmp;
recloc[0] = 0;
recloc[1] = 0;
// std::vector<double> alphap(order);
// std::vector<double> alpham(order);
double sumalpham = 0;
double sumalphap = 0;
if (order == 1){
recloc[0] = uloc[0];
recloc[1] = uloc[0];
}else{
// std::vector<double> umh(order);
for (int r = 0; r < order; r++){
betar = 0;
upl = 0;
uml = 0;
for (int j = 0; j < order; j++){
upl += Crec[r+1][j]*uloc[order-1 - r + j];
uml += Crec[r][j]*uloc[order-1 - r + j];
for (int i = 0; i < order; i++){
betar += uloc[order-1 - r + j]*beta[j][i][r]*uloc[order-1 - r + i];
}
}
tmp = (vareps+betar);
sumalpham += dw[order-1-r]/(tmp*tmp);
recloc[0]+= dw[order-1-r]/(tmp*tmp)*uml;
sumalphap += dw[r]/(tmp*tmp);
recloc[1]+= dw[r]/(tmp*tmp)*upl;
}
recloc[0] = recloc[0]/sumalpham;
recloc[1] = recloc[1]/sumalphap;
// recloc[0] = uloc[1];
// recloc[1] = uloc[1];
}
}
c4science · Help