05.03., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

Commit 8026c92c authored by Anne Reinarz's avatar Anne Reinarz

Cleanup of not-working/maintained applications - moved to branch old_applications

parent 677e1ca8
/**
*
* Euler Flow
*
* This setup corresponds to the Finite Volume demonstrator code as discussed in
* Section 2 of the guidebook. To keep stuff here extremely simple, the spec file
* does not contain any global-optimisations or any parallelisation. How to add parallel
* features is solely described in the guidebook.
*
*
* 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
*/
exahype-project EulerFV
peano-kernel-path const = ./Peano
exahype-path const = ./ExaHyPE
output-directory const = ./ApplicationExamples/CPC/Euler
computational-domain
dimension const = 2
width = 1.0, 1.0
offset = 0.0, 0.0
end-time = 0.1
end computational-domain
solver Finite-Volumes MyEulerSolver
variables const = rho:1,j:3,E:1
patch-size const = 10
maximum-mesh-size = 2e-2
time-stepping = global
type const = godunov
terms const = flux
optimisation const = generic
language const = C
plot vtu::Cartesian::cells::ascii EulerWriter
variables const = 5
time = 0.0
repeat = 0.01
output = ./variables
end plot
end solver
end exahype-project
How to install:
===============
- Build the toolkit (if not done before):
-- Change into Toolkit directory
-- Type
./build.sh
- Generate all gluecode
-- Change into the ExaHyPE root directory
-- Type in
java -jar Toolkit/dist/ExaHyPE.jar Demonstrators/EulerFV/EulerFV.exahype
- Change into Demonstrators/EulerFV
-- If you prefer the GNU compiler, change the default compiler
export COMPILER=gnu
-- Type in
make
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "EulerWriter.h"
EulerFV::EulerWriter::EulerWriter(MyEulerSolver& solver) {
// @TODO Please insert your code here.
}
EulerFV::EulerWriter::~EulerWriter() {
}
void EulerFV::EulerWriter::startPlotting( double time) {
// @TODO Please insert your code here.
}
void EulerFV::EulerWriter::finishPlotting() {
// @TODO Please insert your code here.
}
void EulerFV::EulerWriter::mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp
) {
const int writtenUnknowns = 5;
for (int i=0; i<writtenUnknowns; i++){
outputQuantities[i] = Q[i];
}
}
\ No newline at end of file
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef POSTPROCESSING_EulerWriter_CLASS_HEADER_
#define POSTPROCESSING_EulerWriter_CLASS_HEADER_
#include "exahype/plotters/Plotter.h"
namespace EulerFV {
/* Forward declaration: */
class MyEulerSolver;
class EulerWriter;
}
class EulerFV::EulerWriter : public exahype::plotters::Plotter::UserOnTheFlyPostProcessing {
public:
EulerWriter(MyEulerSolver& solver);
virtual ~EulerWriter();
void startPlotting(double time) override;
void finishPlotting() override;
void mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp) override;
};
#endif /* POSTPROCESSING_EulerWriter_CLASS_HEADER_ */
\ No newline at end of file
This diff is collapsed.
#include "MyEulerSolver.h"
#include "MyEulerSolver_Variables.h"
#include "Logo.h"
tarch::logging::Log EulerFV::MyEulerSolver::_log( "EulerFV::MyEulerSolver" );
void EulerFV::MyEulerSolver::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
// @todo Please implement/augment if required
}
void EulerFV::MyEulerSolver::adjustSolution(const double* const x,const double t,const double dt, double* const Q) {
if ( tarch::la::equals( t,0.0 ) ) {
Variables vars(Q);
tarch::la::Vector<DIMENSIONS,double> myX( x[0] - 0.06, 1.0-x[1] - 0.25 ); // translate
myX *= static_cast<double>(Image.width);
tarch::la::Vector<DIMENSIONS,int> myIntX( 1.2*myX(0) , 1.2*myX(1) ); // scale
double Energy = 0.1;
if (
myIntX(0) > 0 && myIntX(0) < static_cast<int>(Image.width)
&&
myIntX(1) > 0 && myIntX(1) < static_cast<int>(Image.height)
) {
Energy += 1.0-Image.pixel_data[myIntX(1)*Image.width+myIntX(0)];
}
vars.rho() = 1.0;
vars.E() = Energy;
vars.j(0,0,0);
}
}
void EulerFV::MyEulerSolver::eigenvalues(const double* const Q, const int normalNonZeroIndex, double* const lambda) {
ReadOnlyVariables vars(Q);
Variables eigs(lambda);
const double GAMMA = 1.4;
const double irho = 1./vars.rho();
const double p = (GAMMA-1) * (vars.E() - 0.5 * irho * vars.j()*vars.j() );
double u_n = Q[normalNonZeroIndex + 1] * irho;
double c = std::sqrt(GAMMA * p * irho);
eigs.rho()=u_n - c;
eigs.E() =u_n + c;
eigs.j(u_n,u_n,u_n);
}
void EulerFV::MyEulerSolver::flux(const double* const Q, double** const F) {
ReadOnlyVariables vars(Q);
Fluxes fluxes(F);
tarch::la::Matrix<3,3,double> I;
I = 1, 0, 0,
0, 1, 0,
0, 0, 1;
const double GAMMA = 1.4;
const double irho = 1./vars.rho();
const double p = (GAMMA-1) * (vars.E() - 0.5 * irho * vars.j()*vars.j() );
fluxes.rho ( vars.j() );
fluxes.j ( irho * outerDot(vars.j(),vars.j()) + p*I );
fluxes.E ( irho * (vars.E() + p) * vars.j() );
}
void EulerFV::MyEulerSolver::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int normalNonZero,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut) {
ReadOnlyVariables varsInside(stateInside);
Variables varsOutside(stateOutside);
varsOutside = varsInside;
}
#ifndef __MyEulerSolver_CLASS_HEADER__
#define __MyEulerSolver_CLASS_HEADER__
// This file was initially generated by the ExaHyPE toolkit.
// You can modify it in order to extend your solver with features.
// Whenever this file is present, a re-run of the ExaHyPE toolkit will
// not overwrite it. Delete it to get it regenerated.
//
// ========================
// www.exahype.eu
// ========================
#include <ostream>
#include "AbstractMyEulerSolver.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
namespace EulerFV{
class MyEulerSolver;
}
class EulerFV::MyEulerSolver : public EulerFV::AbstractMyEulerSolver {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
MyEulerSolver(const double maximumMeshSize,const exahype::solvers::Solver::TimeStepping timeStepping);
/**
* Initialise the solver.
*
* \param[in] cmdlineargs the command line arguments.
* \param[in] constants access to the constants specified for the solver.
*/
void init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) final override ;
/**
* @see FiniteVolumesSolver
*/
void adjustSolution(const double* const x,const double t,const double dt, double* const Q) override;
/**
* Compute the eigenvalues of the flux tensor per coordinate direction \p d.
*
* \param[in] Q the conserved variables associated with a quadrature node
* as C array (already allocated).
* \param[in] d the column of the flux vector (d=0,1,...,DIMENSIONS).
* \param[inout] lambda the eigenvalues as C array (already allocated).
*/
void eigenvalues(const double* const Q,const int d,double* const lambda) override;
/**
* Impose boundary conditions at a point on a boundary face
* within the time interval [t,t+dt].
*
* \param[in] x the physical coordinate on the face.
* \param[in] t the start of the time interval.
* \param[in] dt the width of the time interval.
* \param[in] faceIndex indexing of the face (0 -- {x[0]=xmin}, 1 -- {x[1]=xmax}, 2 -- {x[1]=ymin}, 3 -- {x[2]=ymax}, and so on,
* where xmin,xmax,ymin,ymax are the bounds of the cell containing point x.
* \param[in] d the coordinate direction the face normal is pointing to.
* \param[in] QIn the conserved variables at point x from inside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
* \param[inout] QOut the conserved variables at point x from outside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
*/
void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int normalNonZero,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut) override;
/**
* Compute the flux tensor.
*
* \param[in] Q the conserved variables (and parameters) associated with a quadrature point
* as C array (already allocated).
* \param[inout] F the fluxes at that point as C array (already allocated).
*/
void flux(const double* const Q,double** const F) override;
/* viscousFlux() function not included, as requested in the specification file */
/* algebraicSource() function not included, as requested by the specification file */
/* nonConservativeProduct() function is not included, as requested in the specification file */
/* pointSource() function not included, as requested in the specification file */
};
#endif // __MyEulerSolver_CLASS_HEADER__
! C2P for GRMHD
! C2P for GRMHD
RECURSIVE SUBROUTINE PDEPrim2Cons(Q,V)
USE Parameters, ONLY: gamma, nDim
IMPLICIT NONE
#ifdef VECTOR
#ifdef AVX512
INTEGER, PARAMETER :: nVar = 24 ! The number of variables of the PDE system
#else
INTEGER, PARAMETER :: nVar = 20 ! The number of variables of the PDE system
#endif
#else
INTEGER, PARAMETER :: nVar = 19 ! The number of variables of the PDE system
#endif
! Argument list declaration
REAL :: Q(nVar), V(nVar)
INTENT(IN) :: V
INTENT(OUT) :: Q
! Local variable declaration
REAL :: rho,p,psi,lapse,gp,gm
REAL :: shift(3),v_cov(3),v_contr(3),vxB_cov(3),vxB_contr(3),B_cov(3),B_contr(3)
REAL :: g_cov(3,3),g_contr(3,3)
REAL :: vb,v2,e2,b2,uem,LF,gamma1,w,ww
INTEGER :: i
rho = V(1)
v_cov = V(2:4)
p = V(5)
!
psi = V(9)
lapse = V(10)
DO i=1,3
B_contr(i) = V(5+i) !wrong
shift(i) = V(10+i) ! NB: we choose V() and Q() being the shift_controvariant!
ENDDO
!
!gammaij = V(14:19)
g_cov(1,1) = V(14)
g_cov(1,2) = V(15)
g_cov(1,3) = V(16)
g_cov(2,2) = V(17)
g_cov(2,3) = V(18)
g_cov(3,3) = V(19)
g_cov(2,1) = V(15)
g_cov(3,1) = V(16)
g_cov(3,2) = V(18)
!
CALL MatrixInverse3x3(g_cov,g_contr,gp)
gp = SQRT(gp)
gm = 1./gp
!
v_contr = MATMUL(g_contr,v_cov)
B_cov = MATMUL(g_cov,B_contr(1:3))
!
! COVARIANT =>
vxB_cov(1) = v_contr(2)*B_contr(3) - v_contr(3)*B_contr(2)
vxB_cov(2) = v_contr(3)*B_contr(1) - v_contr(1)*B_contr(3)
vxB_cov(3) = v_contr(1)*B_contr(2) - v_contr(2)*B_contr(1)
vxB_cov(1) = gp*vxB_cov(1)
vxB_cov(2) = gp*vxB_cov(2)
vxB_cov(3) = gp*vxB_cov(3)
!
! CONTRAVARIANT =>
vxB_contr(1) = v_cov(2)*B_cov(3) - v_cov(3)*B_cov(2)
vxB_contr(2) = v_cov(3)*B_cov(1) - v_cov(1)*B_cov(3)
vxB_contr(3) = v_cov(1)*B_cov(2) - v_cov(2)*B_cov(1)
vxB_contr(1) = gm*vxB_contr(1)
vxB_contr(2) = gm*vxB_contr(2)
vxB_contr(3) = gm*vxB_contr(3)
!
vb = v_cov(1)*B_contr(1) + v_cov(2)*B_contr(2) + v_cov(3)*B_contr(3)
v2 = v_contr(1)*v_cov(1) + v_contr(2)*v_cov(2) + v_contr(3)*v_cov(3)
e2 = vxB_contr(1)*vxB_cov(1) + vxB_contr(2)*vxB_cov(2) + vxB_contr(3)*vxB_cov(3)
b2 = B_contr(1)*B_cov(1) + B_contr(2)*B_cov(2) + B_contr(3)*B_cov(3)
!
uem = 0.5*(b2 + e2)
!
!IF (v2 > 1.0) THEN
! ! IF(MAXVAL(ABS(xGP)).LE.Exc_radius) THEN ! RETURN
! WRITE(*,*)'Superluminal velocity in PDEPrim2Cons!!'
! STOP
!ENDIF
lf = 1.0 / sqrt(1.0 - v2)
gamma1 = gamma/(gamma-1.0)
w = rho + gamma1*p
ww = w*lf**2
!
Q(1) = rho*lf
Q(2) = ww*v_cov(1) + b2*v_cov(1) - vb*B_cov(1)
Q(3) = ww*v_cov(2) + b2*v_cov(2) - vb*B_cov(2)
Q(4) = ww*v_cov(3) + b2*v_cov(3) - vb*B_cov(3)
Q(5) = ww - p + uem - Q(1) !!!!! we subtract PDE(Q(1))!!!!
Q(6) = V(6)
Q(7) = V(7)
Q(8) = V(8)
!
DO i=1,8
Q(i) = gp*Q(i)
ENDDO