Commit 2dd2fd52 authored by Houssein Zaghdane's avatar Houssein Zaghdane

personnel work

parent 144a6e4f
......@@ -9,13 +9,13 @@
"architecture": "hsw",
"computational_domain": {
"dimension": 2,
"end_time": 0.1,
"end_time": 3.0,
"offset": [
0.0,
-10.0,
0.0
],
"width": [
1.0,
70.0,
1.0
]
},
......@@ -32,7 +32,7 @@
"type": "Limiting-ADER-DG",
"name": "MySWESolver",
"order": 1,
"maximum_mesh_size": 0.04,
"maximum_mesh_size": 0.1,
"maximum_mesh_depth": 0,
"time_stepping": "global",
"aderdg_kernel": {
......@@ -86,13 +86,50 @@
],
"plotters": [
{
"type": "vtk::Cartesian::cells::ascii",
"type": "vtu::Cartesian::cells::ascii",
"name": "ConservedWriter",
"time": 0.0,
"repeat": 0.01,
"repeat": 0.1,
"output": "./vtk-output/conserved",
"variables": 5
},
{
"type": "probe::ascii",
"name": "ProbeWriter1",
"time": 0.0,
"repeat": 0.001,
"output": "./buoy/waveheight1",
"variables": 2,
"select": {
"x": 1.5,
"y": 0.5
}
},
{
"type": "probe::ascii",
"name": "ProbeWriter2",
"time": 0.0,
"repeat": 0.001,
"output": "./buoy/waveheight2",
"variables": 2,
"select": {
"x": 2.0,
"y": 0.5
}
},
{
"type": "probe::ascii",
"name": "ProbeWriter3",
"time": 0.0,
"repeat": 0.001,
"output": "./buoy/waveheight3",
"variables": 2,
"select": {
"x": 2.5,
"y": 0.5
}
}
]
}
]
......
......@@ -37,4 +37,13 @@ void SWE::ConservedWriter::mapQuantities(
outputQuantities[i] = Q[i];
}
outputQuantities[4] = Q[0] + Q[3];
for (int i=0; i<writtenUnknowns; i++)
{
if (outputQuantities[i] < 1e-8){
outputQuantities[i] = 0.0;
}
}
}
......@@ -20,14 +20,15 @@
using namespace kernels;
double grav;
double H;
tarch::logging::Log SWE::MySWESolver_ADERDG::_log( "SWE::MySWESolver_ADERDG" );
void SWE::MySWESolver_ADERDG::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
grav = 9.81;
//std::cout << "Param " << muq::param[0] << " " << muq::param[1] << std::endl;
std::cout << "Param " << muq::param[0] << std::endl;
H = muq::param[0];
//for testing csv writer/reader
/*const int nelem = 9*3;
std::vector<double> a((nelem+1)*(nelem+1));
......@@ -84,33 +85,49 @@ void SWE::MySWESolver_ADERDG::adjustSolution(double* const luh, const tarch::la:
kernels::idx2 id_xy(basisSize,basisSize);
int num_nodes = basisSize;
double offset_x=center[0]-0.5*dx[0];
double offset_y=center[1]-0.5*dx[1];
const double d = 0.15;
//H = 0.045;
const double beta = std::atan(1/19.85);
double gamma = std::sqrt((3*H)/(4*d));
double x0 = d * cos(beta)/sin(beta);
double L = d * std::log(std::sqrt(20) + std::sqrt(20 - 1)) / gamma;
for (int i=0; i< num_nodes; i++){
for (int j=0; j< num_nodes; j++){
double x = (offset_x+dx[0]*kernels::gaussLegendreNodes[basisSize-1][i]);
double y = (offset_y+dx[1]*kernels::gaussLegendreNodes[basisSize-1][j]);
double b = 0.0;//SWE::linearInterpolation(x,y,a);
luh[id_xyf(i,j,0)] = initial_waterheight(x,y);
luh[id_xyf(i,j,1)] = 0; //hu
double eta = H * (1/(std::cosh(gamma*(x-(x0 + L))/d))) * (1/(std::cosh(gamma*(x-(x0 + L))/d)));
if (x < 0){
luh[id_xyf(i,j,0)] = 0;
luh[id_xyf(i,j,3)] = -x * std::sin(beta)/std::cos(beta) + d;
}
else if (x <= x0)
{
luh[id_xyf(i,j,0)] = x * std::sin(beta)/std::cos(beta);
luh[id_xyf(i,j,3)] = d - luh[id_xyf(i,j,0)];
}
else{
luh[id_xyf(i,j,0)] = H * (1/(std::cosh(gamma*(x-(x0 + L))/d))) * (1/(std::cosh(gamma*(x-(x0 + L))/d))) + d;
luh[id_xyf(i,j,3)] = 0;
}
luh[id_xyf(i,j,1)] = -eta * std::sqrt(grav/d) * luh[id_xyf(i,j,0)]; //hu
luh[id_xyf(i,j,2)] = 0; //hv
luh[id_xyf(i,j,3)] = b; //b
}
}
}
constexpr int order = MySWESolver_ADERDG::Order;
constexpr int numberOfUnknowns = MySWESolver_ADERDG::NumberOfVariables;
std::vector<std::vector<double>> probe_point = { {0.2,0.2}, {0.4,0.2}, {0.6,0.2}, {0.8,0.2},
{0.2,0.4}, {0.4,0.4}, {0.6,0.4}, {0.8,0.4},
{0.2,0.6}, {0.6,0.6}, {0.8,0.6},
{0.2,0.8}, {0.4,0.8}, {0.6,0.8}, {0.8,0.8},
{0.4,0.6}
std::vector<std::vector<double>> probe_point = { {1.5,0.5}, {2,0.5}, {2.5,0.5} //, {3, 0.5}, {4, 0.5}
};
for (int i = 0; i< probe_point.size(); i++){
if( isInside(probe_point[i], center, dx)){
......@@ -118,6 +135,7 @@ void SWE::MySWESolver_ADERDG::adjustSolution(double* const luh, const tarch::la:
double center_[DIMENSIONS] ={center[0],center[1]};
double dx_[DIMENSIONS] ={dx[0],dx[1]};
muq::solution[i] = interpolate( center_, dx_, &(probe_point[i][0]), numberOfUnknowns, 0, order, luh);
std::cout << "solution in aderdg " << muq::solution[i] << std::endl;
}
}
}
......@@ -157,7 +175,7 @@ void SWE::MySWESolver_ADERDG::eigenvalues(const double* const Q,const int d,doub
// Number of variables + parameters = 4 + 0
ReadOnlyVariables vars(Q);
Variables eigs(lambda);
const double c = std::sqrt(grav*vars.h());
const double ih = 1./vars.h();
double u_n = Q[d + 1] * ih;
......@@ -204,7 +222,7 @@ void SWE::MySWESolver_ADERDG::flux(const double* const Q,double** const F) {
}
else {
f[0] = vars.hu();
f[1] = vars.hu() * vars.hu() * ih + 0.5 * grav * vars.h() * vars.h();
f[1] = vars.hu() * vars.hu() * ih + + 0.5 * grav * vars.h() * vars.h();;
f[2] = vars.hu() * vars.hv() * ih;
f[3] = 0.0;
......
......@@ -12,6 +12,8 @@ int scenario_FV;
tarch::logging::Log SWE::MySWESolver_FV::_log( "SWE::MySWESolver_FV" );
void SWE::MySWESolver_FV::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
epsilon_FV = 10.00 / 1000.0;
grav_FV = 9.81 / 1000.0;
}
......@@ -20,7 +22,34 @@ void SWE::MySWESolver_FV::adjustSolution(const double* const x,const double t,co
// Number of variables = 4 + #parameters
if (tarch::la::equals(t,0.0)) {
std::cout << "Warning: using FV in initial step" << std::endl;
//std::cout << "Warning: using FV in initial step" << std::endl;
/*MySWESolver_FV::Variables vars(Q);
const double d = 0.15;
const double H = 0.3 * d;
const double beta = std::atan(1/19.85);
double gamma = std::sqrt((3*H)/(4*d));
double x0 = d * cos(beta)/sin(beta);
double L = d * std::log(std::sqrt(20) + std::sqrt(20 - 1)) / gamma;
double eta = H * (1/(std::cosh(gamma*(x[0]-(x0 + L))/d))) * (1/(std::cosh(gamma*(x[0]-(x0 + L))/d)));
if (x[0] < 0){
vars.h() = 0;
vars.b() = -x[0] * std::sin(beta)/std::cos(beta) + d;
}
else if (x[0] <= x0){
vars.h() = x[0] * std::sin(beta)/std::cos(beta);
vars.b() = d - vars.h();
}
else{
vars.h() = H * (1/(std::cosh(gamma*(x[0]-(x0 + L))/d))) * (1/(std::cosh(gamma*(x[0]-(x0 + L))/d))) + d;
vars.b() = 0;
}
vars.hu() = -eta * std::sqrt(grav_FV/d) * vars.h();
vars.hv() = 0.0;*/
}
else{
if(Q[0] < epsilon_FV){
......@@ -33,7 +62,6 @@ void SWE::MySWESolver_FV::adjustSolution(const double* const x,const double t,co
void SWE::MySWESolver_FV::eigenvalues(const double* const Q, const int dIndex, double* const lambda) {
// Dimensions = 2
// Number of variables = 4 + #parameters
ReadOnlyVariables vars(Q);
Variables eigs(lambda);
......@@ -44,6 +72,7 @@ void SWE::MySWESolver_FV::eigenvalues(const double* const Q, const int dIndex, d
eigs.hu() = u_n - c;
eigs.hv() = u_n;
eigs.b() = 0.0;
}
void SWE::MySWESolver_FV::boundaryValues(
......
......@@ -247,11 +247,10 @@ int init(int argc, char** argv) {
//2. registerSolvers(), do I need to run this again?
std::vector<double> run_exahype(std::vector<double> param_, int level){
muq::param.resize(2);
muq::solution.resize(16);
muq::param.resize(1);
muq::solution.resize(3);
param=param_; //store parameters
mesh_width = std::pow(3,-level);
std::cout << "Mesh size " << mesh_width << std::endl;
//mesh_width = std::pow(3,-level);
if(!firstRun){
tarch::parallel::NodePool::getInstance().shutdown();
tarch::parallel::NodePool::getInstance().init();
......@@ -263,6 +262,7 @@ int init(int argc, char** argv) {
firstRun = false;
runner.shutdownHeaps();
peano::parallel::JoinDataBufferPool::getInstance().releaseMessages();
return solution;
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment