In January 2021 we will introduce a 10 GB quota for project repositories. Higher limits for individual projects will be available on request. Please see https://doku.lrz.de/display/PUBLIC/GitLab for more information.

Commit fe55a9b7 authored by leo.rannabauer's avatar leo.rannabauer

Init code frame fro geoclaw SWE solver

parent 71bbd773
/**
* This file is part of the ExaHyPE project.
* Copyright (c) 2016 http://exahype.eu
* All rights reserved.
*
* The project has received funding from the European Union's Horizon
* 2020 research and innovation programme under grant agreement
* No 671698. For copyrights and licensing, please consult the webpage.
*
* Released under the BSD 3 Open Source License.
* For the full license text, see LICENSE.txt
**/
/**
2D-SWE with Godunov FV solver.
*/
exahype-project SWE
peano-kernel-path const = ./Peano
exahype-path const = ./ExaHyPE
output-directory const = ./ApplicationExamples/SWE/SWE_geoclaw
architecture const = noarch
//log-file = mylogfile.log
computational-domain
dimension const = 2
/* for run-up scenarios in channels */
/* width = 120.0, 1.0 */
/* for shelf scenarios */
/*width = 1200.0, 1200.0 */
/* for dam breaks */
width = 4.0, 4.0
offset = -2.0, -2.0
end-time = 5.0
end computational-domain
shared-memory
identifier = dummy
configure = {background-tasks:8}
cores = 8
properties-file = sharedmemory.properties
end shared-memory
global-optimisation
fuse-algorithmic-steps = all
fuse-algorithmic-steps-rerun-factor = 0.99
fuse-algorithmic-steps-diffusion-factor = 0.99
spawn-predictor-as-background-thread = off
spawn-amr-background-threads = off
/* 0.0 und 0.8 sind schon mal zwei Faktoren */
disable-vertex-exchange-in-time-steps = on
time-step-batch-factor = 0.0
disable-metadata-exchange-in-batched-time-steps = off
double-compression = 0.0
spawn-double-compression-as-background-thread = off
end global-optimisation
solver Finite-Volumes MySWESolver
variables const = h:1,hu:1,hv:1,b:1
//parameters const = b:1
patch-size const = 3
maximum-mesh-size = 0.05
time-stepping = global
type const = godunov
terms const = flux
optimisation const = generic, usestack
language const = C
constants = grav:9.81, epsilon:1e-3, scenario:11
plot vtu::Cartesian::cells::ascii ConservedWriter
variables const = 6
time = 0.0
repeat = 0.1
output = ./output/conserved
end plot
end solver
end exahype-project
/*------Scenarios-------
0: ShockShockProblem (width = 10.0, 10.0 offset = 0.0, 0.0)
1: RareRareProblem (width = 10.0, 10.0 offset = 0.0, 0.0)
2: GaussFunction (width = 10.0, 10.0 offset = 0.0, 0.0)
3: ExpBreakProblem
4: DamBreakProblem (width = 10.0, 10.0 offset = 0.0, 0.0)
5: SeaAtRestProblem
6: SteadyRunUpLinear
7: RunUpLinear
8: SteadyRunUpShelf
9: RunUpShelf
10: WettingDryingProblem (width = 20.0, 1.0 offset = -10.0, 0.0)
11: OscillatingLake (width = 4.0, 4.0 offset = -2.0, -2.0)
12: RunUpTest (width = 10.0, 10.0 offset = 0.0, 0.0)
13: SolitaryWaveOnSimpleBeach (Case = 0.0185) (width = 70.0, 1.0 offset = -10.0, 0.0)
14: SolitaryWaveOnSimpleBeach (Case = 0.3) (width = 70.0, 1.0 offset = -10.0, 0.0)
Default: GaussFunction (width = 10.0, 10.0 offset = 0.0, 0.0)*/
#include "InitialData.h"
#include "MySWESolver_Variables.h"
#include "MySWESolver.h"
#include <cmath>
using namespace std;
///// 2D /////
extern double grav;
extern int scenario;
#ifdef Dim2
/*
* Constant water height with both sides smashing together
*/
void SWE::ShockShockProblem(const double * const x, double* Q){
MySWESolver::Variables vars(Q);
if(x[0] < 5) {
vars.h() = 4.0;
vars.hu()= 2.0;
vars.hv()= 0.0;
vars.b() = 0;
} else {
vars.h() = 4.0;
vars.hu()= -2.0;
vars.hv()= 0.0;
vars.b() = 0.0;
}
}
/*
* Constant water height with both sides moving away from each other
*/
void SWE::RareRareProblem(const double * const x, double* Q){
MySWESolver::Variables vars(Q);
if(x[0] < 5) {
vars.h() = 4.0;
vars.hu()= -2.0;
vars.hv()= 0.0;
vars.b() = 0;
} else {
vars.h() = 4.0;
vars.hu()= 2.0;
vars.hv()= 0.0;
vars.b() = 0.0;
}
}
/*
* Gaussfunktion
*/
void SWE::GaussFunctionProblem(const double* const x, double* Q){
MySWESolver::Variables vars(Q);
vars.h() = exp(-pow(x[0] - 5, 2)) + 1;
vars.hu() = 0.0;
vars.hv() = 0.0;
vars.b() = 0.0;
}
/*
* Constant water height with "exponential" hump in bathymetry.
*/
void SWE::ExpBreakProblem(const double* const x,double* Q) {
MySWESolver::Variables vars(Q);
if((x[0] -5) *(x[0] -5) + (x[1] -5) *(x[1] -5) < 2) {
vars.h() = 4.0;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = exp(-pow((x[0]-5),2)-pow((x[1]-5),2));
} else {
vars.h() = 4.0;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = 0.0;
}
}
/*
* Constant water height with discontinuous jump in bathymetry.
*/
void SWE::DamBreakProblem(const double* const x,double* Q) {
MySWESolver::Variables vars(Q);
if((x[0] -5) *(x[0] -5) + (x[1] -5) *(x[1] -5) < 2) {
vars.h() = 4.0;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = 1;
} else {
vars.h() = 4.0;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = 0.0;
}
}
/*
* Sea at rest (steady state).
*/
void SWE::SeaAtRestProblem(const double* const x,double* Q) {
MySWESolver::Variables vars(Q);
if(x[0] >= (10.0/3.0) && x[0]<= (20.0/3.0)) {
vars.h() = 2 - x[0]*(3.0/10.0) + 1;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = x[0]*(3.0/10.0) - 1;
} else {
vars.h() = 2.0;
vars.hu()= 0.0;
vars.hv()= 0.0;
vars.b() = 0.0;
}
}
/*
* Simulates a channel with a linearly rising ramp at the left end of the channel.
* To be used with WALL boundaries.
*/
void SWE::SteadyRunUpLinear(const double* const x, double* Q) {
MySWESolver::Variables vars(Q);
const double d =1;
const double xr=19.85;
vars.b() = (x[0]>xr) ? 0 : (-d/xr)*x[0] +d;
vars.h() = (x[0]>xr) ? d : (d/xr)*x[0];
vars.hu() = 0.0;
vars.hv() = 0.0;
}
/*
* Simulates a wave which approaches a linearly rising ramp at the left end of the channel.
* To be used with WALL boundaries.
*/
void SWE::RunUpLinear(const double* const x, double* Q) {
MySWESolver::Variables vars(Q);
const double as=0.3;
const double d =1;
const double xr=19.85;
const double alpha=atan(d/xr);
const double xs= d/tan(alpha) + sqrt(4.0/(3.0*as))*acosh(sqrt(20.0));
vars.b() = (x[0]>xr) ? 0 : (-d/xr)*x[0] +d;
vars.h() = (x[0]>xr) ? d : (d/xr)*x[0];
vars.h() += as*(1.0/pow(cosh(sqrt(3.0*as/4.0)*(x[0]-xs)),2));
vars.hu() = -vars.h()* as*(1.0/pow(cosh(sqrt(3*as/4)*(x[0]-xs)),2))*sqrt(grav/d);
vars.hv() = 0.0;
}
/*
* Steady-state test case for artificial continental shelf.
*/
void SWE::SteadyRunUpShelf(const double* const x, double* Q) {
MySWESolver::Variables vars(Q);
const double d1 = 200.0;
const double d2 = 2000.0;
const double xb = 10.0;
const double xe = 50.0;
if(x[0]>xe)
{
vars.h() = d2;
vars.b() = 0.0;
}
else if ( x[0]>xb)
{
vars.h() = ((d1-d2)/(xb-xe))*(x[0]-xb) +d1;
vars.b() = d2 - vars.h();
}
else
{
vars.h() = d1;
vars.b() = d2-d1;
}
vars.hu() = 0.0;
vars.hv() = 0.0;
}
/*
* Run-up on artificial continental shelf.
*/
void SWE::RunUpShelf(const double* const x, double* Q) {
MySWESolver::Variables vars(Q);
const double d1 = 200.0;
const double d2 = 2000.0;
const double xb = 10.0;
const double xe = 50.0;
const double xw = 500.0;
const double aw = 1.0;
vars.hu()=0.0;
vars.hv()=0.0;
if(x[0]>xe)
{
vars.h() = d2;
vars.b() = 0.0;
}
else if (x[0]>xb)
{
vars.h() = ((d1-d2)/(xb-xe))*(x[0]-xb) +d1;
vars.b() = d2 - vars.h();
}
else
{
vars.h() = d1;
vars.b() = d2-d1;
}
vars.h() += aw*(1.0/pow(cosh(x[0]-xw),2));
vars.hu() +=-vars.h()*aw*(1.0/pow(cosh(x[0]-xw),2))*sqrt(grav/d2);
//vars.hu() = -vars.h()* aw*(1.0/pow(cosh(sqrt(x[0]-xw),2))*sqrt(grav/d);
//vars.hu() = -vars.h()* sqrt(grav*vars.h());
}
// width = 10.0,10.0
// offset = 0.0, 0.0
void SWE::WettingDryingProblem(const double* const x, double* Q){
MySWESolver::Variables vars(Q);
if(x[0] < -5) {
vars.h() = 2.0;
} else {
vars.h() = 0.0;
}
vars.hu() = 0.0;
vars.hv() = 0.0;
vars.b() = -0.1*x[0] + exp((-x[0]*x[0])/(0.1*0.1)) + 20;
}
void SWE::OscillatingLake(const double* const x, double* Q){
MySWESolver::Variables vars(Q);
double omega = sqrt(0.2*grav);
vars.b() = 0.1 * (pow(x[0], 2) + pow(x[1], 2));
vars.h() = max(0.0, 0.05 * (2 * x[0] * cos(omega * 0) + 2 * x[1] * sin(omega * 0)) + 0.075 - vars.b());
vars.hu() = 0.5 * omega * sin(omega * 0) * vars.h();
vars.hv() = 0.5 * omega * cos(omega * 0) * vars.h();
}
// width = 10.0,10.0
// offset = 0.0, 0.0
void SWE::RunUpTest(const double* const x, double* Q){
MySWESolver::Variables vars(Q);
if(x[0] < 7){
vars.h() = 0.5*exp(-pow(x[0] - 3.5, 2)) + 0.5;
vars.hu() = 2.0;
vars.b() = 0.0;
}
else {
vars.h() = 0.0;
vars.hu() = 0.0;
vars.b() = (x[0] - 7) * 1;
}
vars.hv() = 0.0;
}
//width = 70.0,10.0
//offset = -10.0,0.0
void SWE::SolitaryWaveOnSimpleBeach(const double*const x, double* Q, double Hd, double d){
MySWESolver::Variables vars(Q);
const double H = Hd * 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/d) * vars.h();
vars.hv() = 0.0;
}
#endif
void SWE::initialData(const double* const x,double* Q) {
switch (scenario)
{
case 0:
ShockShockProblem(x, Q);
break;
case 1:
RareRareProblem(x, Q);
break;
case 2:
GaussFunctionProblem(x, Q);
break;
case 3:
ExpBreakProblem(x, Q);
break;
case 4:
DamBreakProblem(x, Q);
break;
case 5:
SeaAtRestProblem(x, Q);
break;
case 6:
SteadyRunUpLinear(x, Q);
break;
case 7:
RunUpLinear(x, Q);
break;
case 8:
SteadyRunUpShelf(x, Q);
break;
case 9:
RunUpShelf(x, Q);
break;
case 10:
WettingDryingProblem(x, Q);
break;
case 11:
OscillatingLake(x, Q);
break;
case 12:
RunUpTest(x, Q);
break;
case 13:
SolitaryWaveOnSimpleBeach(x, Q, 0.0185, 0.3);
break;
case 14:
SolitaryWaveOnSimpleBeach(x, Q, 0.3, 0.15);
break;
default:
GaussFunctionProblem(x, Q);
}
}
#ifndef __InitialData_CLASS_HEADER__
#define __InitialData_CLASS_HEADER__
namespace SWE {
void ShockShockProblem(const double* const x,double* Q);
void RareRareProblem(const double* const x,double* Q);
void GaussFunctionProblem(const double* const x,double* Q);
void ExpBreakProblem(const double* const x,double* Q);
void DamBreakProblem(const double* const x,double* Q);
void SeaAtRestProblem(const double* const x,double* Q);
void SteadyRunUpLinear(const double* const x,double* Q);
void RunUpLinear(const double* const x,double* Q);
void SteadyRunUpShelf(const double* const x,double* Q);
void RunUpShelf(const double* const x,double* Q);
void WettingDryingProblem(const double* const x, double* Q);
void RunUpTest(const double* const x, double* Q);
void OscillatingLake(const double* const x, double* Q);
void SolitaryWaveOnSimpleBeach(const double*const x, double* Q, double Hd, double d);
void initialData(const double* const x,double* Q);
}
#endif // __InitialData_CLASS_HEADER__
#include "MySWESolver.h"
#include "InitialData.h"
#include "MySWESolver_Variables.h"
#include "kernels/KernelUtils.h"
using namespace kernels;
double grav;
double epsilon;
int scenario;
tarch::logging::Log SWE::MySWESolver::_log( "SWE::MySWESolver" );
void SWE::MySWESolver::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
if (constants.isValueValidDouble( "grav" )) {
grav = constants.getValueAsDouble("grav");
}
if (constants.isValueValidDouble( "epsilon" )) {
epsilon = constants.getValueAsDouble( "epsilon" );
}
if (constants.isValueValidInt( "scenario" )) {
scenario = constants.getValueAsInt( "scenario" );
}
}
void SWE::MySWESolver::adjustSolution(const double* const x,const double t,const double dt, double* const Q) {
// Dimensions = 2
// Number of variables = 4 + #parameters
if (tarch::la::equals(t,0.0)) {
initialData(x, Q);
}
else{
if(Q[0] < epsilon){
Q[1] = 0;
Q[2] = 0;
}
}
}
void SWE::MySWESolver::eigenvalues(const double* const Q, const int dIndex, double* const lambda) {
// Dimensions = 2
// Number of variables = 4 + #parameters
ReadOnlyVariables vars(Q);
Variables eigs(lambda);
const double c = std::sqrt(grav * vars.h());
double u_n = Q[dIndex + 1] * vars.h() * std::sqrt(2)/std::sqrt(std::pow(vars.h(), 4) + std::pow(std::max(vars.h(), epsilon), 4));
eigs.h() = u_n + c;
eigs.hu() = u_n - c;
eigs.hv() = u_n;
eigs.b() = 0.0;
}
void SWE::MySWESolver::boundaryValues(
const double* const x,
const double t,const double dt,
const int faceIndex,
const int d,
const double* const stateInside,
double* const stateOutside) {
// Dimensions = 2
// Number of variables = 4 + #parameters
//OL Code
// stateOutside[0] = stateInside[0];
// stateOutside[1] = 0.0;
// stateOutside[2] = 0.0;
// stateOutside[3] = 0.0;
//normal Code
stateOutside[0] = stateInside[0];
stateOutside[1] = stateInside[1];
stateOutside[2] = stateInside[2];
stateOutside[3] = stateInside[3];
//Wall
stateOutside[d + 1] = -stateInside[d + 1];
}
//***********************************************************
//*********************** PDE *******************************
//***********************************************************
//to add new PDEs specify them in the specification file, delete this file and its header and rerun the toolkit
void SWE::MySWESolver::flux(const double* const Q,double** const F) {
// Dimensions = 2
// Number of variables + parameters = 4 + 0
ReadOnlyVariables vars(Q);
double* f = F[0];
double* g = F[1];
double u_n = vars.hu() * vars.h() *std::sqrt(2)/std::sqrt(std::pow(vars.h(), 4) + std::pow(std::max(vars.h(), epsilon), 4));
double v_n = vars.hv() * vars.h() *std::sqrt(2)/std::sqrt(std::pow(vars.h(), 4) + std::pow(std::max(vars.h(), epsilon), 4));
f[0] = vars.h() * u_n;
f[1] = vars.h() * u_n * u_n; // 0.5 * grav * vars.h() * vars.h();
f[2] = vars.h() * u_n * v_n;
f[3] = 0.0;
g[0] = vars.h() * v_n;
g[1] = vars.h() * u_n * v_n;
g[2] = vars.h() * v_n * v_n; // 0.5 * grav * vars.h() * vars.h();
g[3] = 0.0;
}
// double SWE::MySWESolver::riemannSolver(double* fL, double *fR, const double* qL, const double* qR, const double* gradQL, const double* gradQR, const double* cellSize, int direction) {
// double LL[NumberOfVariables] = {0.0};
// double LR[NumberOfVariables] = {0.0};
// eigenvalues(qL, direction, LL);
// eigenvalues(qR, direction, LR);
// double smax = 0.0;
// for (int i = 0; i < NumberOfVariables; i++) {
// const double abs_sL_i = std::abs(LL[i]);
// smax = std::max( abs_sL_i, smax );
// }
// for (int i = 0; i < NumberOfVariables; i++) {
// const double abs_sR_i = std::abs(LR[i]);
// smax = std::max( abs_sR_i, smax );
// }