Homec4science
Diffusion ConSol (master) Browse

ConSol/ConSol/ConSol/Mesh/BC/SlipWallBC.cpp
f9c1f927e1e5master

Last Change
Raw File

SlipWallBC.cpp
View Options

//
// Created by Fabian Moenkeberg on 2019-05-01.
//
#include "SlipWallBC.hpp"
//#include "flux_farfield.cpp"
void SlipWallBC::initialize(MeshBase *mesh, Configuration config) {
gamma = 1.4;
U_in.resize(4);
U_out.resize(4);
double r_in = 1;
double u1_in = std::sqrt(gamma);
double u2_in =0;
double p_in = 1;
double p_out = 1.3;
double r_out = r_in*(gamma-1+(gamma+1)*p_out)/(gamma+1+(gamma-1)*p_out);
double u1_out = std::sqrt(gamma) + std::sqrt(2)*(1-p_out)/std::sqrt(gamma-1+p_out*(gamma+1));
double u2_out = 0;
double E_in = p_in/(gamma - 1) + 0.5*(u1_in*u1_in+u2_in*u2_in)*r_in;
double E_out = p_out/(gamma - 1) + 0.5*(u1_out*u1_out+u2_out*u2_out)*r_out;
U_in = {r_in, r_in*u1_in, r_in*u2_in, E_in};
U_out = {r_out, r_out*u1_out, r_out*u2_out, E_out};
// U_in = U_out;
// U_out = U_in;
}
void SlipWallBC::updateBoundary(std::vector<std::vector<double> > &U, MeshBase *mesh, double t){
// std::cout << "testB0\n";
// for (int k = 0; k < sysDim; k++){
// std::cout << "testB"<< k<<"\n";
// BC[k].updateBoundary(U, mesh);
// }
// int nCellsTot = mesh->getNcellsTot();
//// std::cout << "testB1\n";
// for (int i = 0; i < nCellsTot; i++){
//// std::cout << "testBi\n";
// U[i] = U[mesh->getBCMap(i)];
// }
std::vector<double> n;
int nBC = mesh->getNbc();
for (int i = 0; i < nBC; i++){
U[mesh->getBCMap(3*i+0)] = U[mesh->getBCMap(3*i+1)];
int iEdge = mesh->getBCMap(3*i+2);
if (iEdge == -1){ // internal hybrid update
U[mesh->getBCMap(3 * i + 0)] = U[mesh->getBCMap(3 * i + 1)];
}else {
std::vector<int> edgesAtBC = mesh->getEdgesAtBC(iEdge);
switch (edgesAtBC[1]) {
case 2: { // slipwall BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
flux_slipWall(U[mesh->getBCMap(3 * i + 1)], n, U[mesh->getBCMap(3 * i + 0)]);
break;
}
case 0: { // means inflow BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
U[mesh->getBCMap(3 * i + 0)] = flux_farfield(U[mesh->getBCMap(3 * i + 1)], n, U_in);
break;
}
case 1: { // means outflow BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
U[mesh->getBCMap(3 * i + 0)] = flux_farfield(U[mesh->getBCMap(3 * i + 1)], n, U_out);
break;
}
}
}
}
int nBC2 = mesh->getNbc2()/2;
int dimSys = mesh->getNelem();
std::vector<double> tmp;
std::vector<double> tmp2;
for (int i = 0; i < nBC2; i++){
tmp = U[mesh->getBCMap2(3*i + 1)];
tmp2 = U[mesh->getBCMap2(3*(i + nBC2) + 1)];
for (int j = 0; j < dimSys; j++){
U[mesh->getBCMap2(3*i + 0)][j] = 0.5*(tmp[j] + tmp2[j]);
}
}
}
void SlipWallBC::updateEdges(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, MeshBase *mesh, double t){
std::vector<std::vector<int>> edgesAtBC = mesh->getEdgesAtBC();
int nelem = mesh->getNelem();
int nQuad = UReconstr[0][0].size();
int nEdgesBC = edgesAtBC.size();
std::vector<double> n;
for(int i = 0; i < nEdgesBC; i++) {
nQuad = UReconstr[0][edgesAtBC[i][0]].size();
for (int k = 0; k < nQuad; k++){
// if ((edgesAtBC[i][1] == 1)||(edgesAtBC[i][1] == 3)){// top and bottom wall with slipwall BC
if ((edgesAtBC[i][1] == 2)){// slipwall BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
flux_slipWall(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
flux_slipWall(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
}
}else if ((edgesAtBC[i][1] == 0)){ // means inflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_in);
}
// }else{ //right wall, means outflow BC
}else if ((edgesAtBC[i][1] == 1)){ // means outflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_out);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_out);
}
}
}
if ((edgesAtBC[i][1] == 11)){ // means hybrid Connection BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
for (int k = 0; k < nQuad; k++){
UReconstr[0][edgesAtBC[i][0]][k] = UReconstr[0][edgesAtBC[i][3]][0];
}
}else{//means left side of edge is out & normal is INWARDS pointing
for ( int l = 0; l < nelem; l++){
UReconstr[1][edgesAtBC[i][0]][0][l] = 0;
for (int k = 0; k < nQuad; k++){
UReconstr[1][edgesAtBC[i][0]][0][l] += UReconstr[1][edgesAtBC[i][3]][k][l]/nQuad;
}
}
}
}else if ((edgesAtBC[i][1] == 14)){ // means hybrid QUADQUAD 2 QUAD or QUAD 2 QUADQUAD Connection BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
UReconstr[0][edgesAtBC[i][0]]= UReconstr[1][edgesAtBC[i][3]];
}else{//means left side of edge is out & normal is INWARDS pointing
UReconstr[1][edgesAtBC[i][0]]= UReconstr[0][edgesAtBC[i][3]];
}
}
}
}
void SlipWallBC::updateEdges(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, MeshBase *mesh, double t, int i){
std::vector<std::vector<int>> edgesAtBC = mesh->getEdgesAtBC();
int nQuad = UReconstr[0][0].size();
int nEdgesBC = edgesAtBC.size();
std::vector<double> n;
for (int k = 0; k < nQuad; k++){
if ((edgesAtBC[i][1] == 2)){// slipwall BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
flux_slipWall(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
flux_slipWall(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
}
}else if ((edgesAtBC[i][1] == 0)){ // means inflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_in);
}
}else if ((edgesAtBC[i][1] == 1)){ // means outflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_out);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_out);
}
}
}
}
void SlipWallBC::flux_slipWall(std::vector<double> U, std::vector<double> n, std::vector<double> &ret){
ret = U;
std::vector<double> V_in{U[1]/U[0],U[2]/U[0]};
double vn = V_in[0]*n[0] + V_in[1]*n[1];
ret[1] = U[0]*(V_in[0]-2*vn*n[0]);
ret[2] = U[0]*(V_in[1]-2*vn*n[1]);
}
c4science · Help