Homec4science
Diffusion ConSol (master) Browse

ConSol/ConSol/ConSol/Reconstr/HybridRBFWENO.cpp
525ec6b60706master

Last Change
Raw File

HybridRBFWENO.cpp
View Options

//
// Created by Fabian Moenkeberg on 2019-10-21.
//
#include "HybridRBFWENO.hpp"
#include "ENO_2D_otf.hpp"
#include "MRWENO_c.hpp"
#include "CWENO_2D.h"
#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 HybridRBFWENO::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Source* source, Configuration config){
if (config.getStencilChoice() == SYMMc){
this->ENO = new class MRWENO_c();
}else if (config.getStencilChoice() == SYMM){
this->ENO = new class CWENO_2D();
}else{
this->ENO = new class ENO_2D_otf();
}
ENO->initialize(model, interpBase, limiter, mesh, source, config);
this->model = model;
this->source = source;
this->limiter = limiter;
this->interpBase = interpBase;
this->order = ceil((config.getOrder()+1)/2.0)+1;
this->nThreads = config.getNthreads();
int n = config.getOrder();
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 (config.getOrder() < 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";
xCoord = mesh->getXCoord();
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];
}
}
}
nCells = mesh->getNcells();
nQuads = QWeights.size();
nDomainsQUAD = mesh->getNdomainsQUAD() + mesh->getNdomainsQUADQUAD() + mesh->getNdomainsRQ();
// nDomainsQUADQUAD = mesh->getNdomainsQUADQUAD();
// nDomainsRQ = mesh->getNdomainsRQ();
nDomainsTRIQUAD = mesh->getNdomainsTRIQUAD();
nDomainsTRI = mesh->getNdomainsTRI();
nDomainsRT = mesh->getNdomainsRT();
submeshCell0 = mesh->getSubmeshCell0();
submeshEdg0 = mesh->getSubmeshEdg0();
nX = mesh->getNx();
nY = mesh->getNy();
// Quadrilateral weights
order_scheme = config.getOrder();
w0.resize(n/2 + (n%2));
abs.resize(n/2 + (n%2));
x_legendre.resize(n);
weights_legendre.resize(n);
w = new double[n];
x = new double[n];
legendre_set(n,x,w);
if (order_scheme < 2){
x_legendre[0] = 0;
weights_legendre[0] = 2;
}else {
for (int i = 0; i < n; i++) {
x_legendre[i] = x[i];
weights_legendre[i] = w[i];
}
}
delete [] w;
delete [] x;
}
void HybridRBFWENO::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U,
MeshBase *mesh, double t, std::vector<std::vector<double>> &Source){
int dimSyst = mesh->getNelem();
ENO->reconstruct(UReconstr, U, mesh, t, Source);
for (int ii = 0; ii < nDomainsQUAD; ii++){
omp_set_num_threads(nThreads);
#pragma omp parallel for
for (int i = 0; i < nX[ii]; i++) {
std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int j = 0; j < nY[ii]-1; j++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 1 - order; l < order; l++){
uloc[l + order - 1] = U[i+ngc+(l+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
}
WENO_ev(uloc, rec);
uRec[k] = rec;
// UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
// UReconstr[1][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
int j = nY[ii]-1;
for (int k = 0; k < dimSyst; k++) {
for (int l = 1 - order; l < order; l++){
uloc[l + order - 1] = U[i+ngc+(l+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
}
WENO_ev(uloc, rec);
uRec[k] = rec;
// UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
// UReconstr[0][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
omp_set_num_threads(nThreads);
#pragma omp parallel for
for (int j = 0; j < nY[ii]; j++) {
std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
std::vector<double> uloc(2*(order-1) + 1);
std::vector<double> rec(2);
for (int i = 0; i < nX[ii]-1; i++) {
for (int k = 0; k < dimSyst; k++) {
for (int l = 1 - order; l < order; l++){
uloc[l+order-1] = U[i + ngc + l + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
}
WENO_ev(uloc, rec);
uRec[k] = rec;
// UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
// UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
int i = nX[ii] - 1;
for (int k = 0; k < dimSyst; k++) {
for (int l = 1 - order; l < order; l++){
uloc[l+order-1] = U[i + ngc + l + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
}
WENO_ev(uloc, rec);
uRec[k] = rec;
// UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
// UReconstr[0][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
}
limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
for (int k = 0; k < dimSyst; k++) {
UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
}
}
}
mesh->updateEdges(UReconstr, t);
}
std::vector<int> HybridRBFWENO::getStencil(int cellIndex, MeshBase *mesh, std::vector<std::vector<double>> *U, int k){
int indMinSmInd = 0;
std::vector<std::vector<int>> stencils = mesh->getStencils(cellIndex);
return stencils[indMinSmInd];
}
void HybridRBFWENO::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];
}
}
void HybridRBFWENO::exportStencils(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, std::string output_file){
ENO->exportStencils(U, mesh, t, output_file);
}
c4science · Help