...
 
Commits (33)
......@@ -11,6 +11,7 @@
[submodule "Submodules/libxsmm"]
path = Submodules/libxsmm
url = https://github.com/hfp/libxsmm.git
branch = release
ignore = dirty #deleted documentation and sample + build
[submodule "Submodules/attrs"]
path = Submodules/attrs
......
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "ConservedWriter.h"
ElasticWave::ConservedWriter::ConservedWriter(ElasticWave::MyElasticWaveSolver& solver) {
// @TODO Please insert your code here.
}
ElasticWave::ConservedWriter::~ConservedWriter() {
}
void ElasticWave::ConservedWriter::startPlotting( double time) {
// @TODO Please insert your code here.
}
void ElasticWave::ConservedWriter::finishPlotting() {
// @TODO Please insert your code here.
}
void ElasticWave::ConservedWriter::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* const Q,
double* const outputQuantities,
double timeStamp
) {
const int writtenUnknowns = 7;
for (int i=0; i<writtenUnknowns; i++){
outputQuantities[i] = Q[i];
}
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef POSTPROCESSING_ConservedWriter_CLASS_HEADER_
#define POSTPROCESSING_ConservedWriter_CLASS_HEADER_
#include "exahype/plotters/Plotter.h"
namespace ElasticWave {
class MyElasticWaveSolver;
class ConservedWriter;
}
class ElasticWave::ConservedWriter : public exahype::plotters::Plotter::UserOnTheFlyPostProcessing {
public:
ConservedWriter(ElasticWave::MyElasticWaveSolver& solver);
virtual ~ConservedWriter();
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* const Q,
double* const outputQuantities,
double timeStamp) override;
};
#endif /* POSTPROCESSING_ConservedWriter_CLASS_HEADER_ */
\ No newline at end of file
{
"project_name": "ElasticWave",
"paths": {
"peano_kernel_path": "./Peano",
"exahype_path": "./ExaHyPE",
"output_directory": "./ApplicationExamples/AdjointTest2D"
},
"architecture": "noarch",
"computational_domain": {
"dimension": 2,
"end_time": 1.0,
"offset": [
0.0,
0.0
],
"width": [
10.0,
10.0
]
},
"solvers": [
{
"type": "ADER-DG",
"name": "MyElasticWaveSolver",
"order": 3,
"maximum_mesh_size": 0.5,
"time_stepping": "globalfixed",
"aderdg_kernel": {
"language": "C",
"nonlinear": false,
"terms": [
"flux",
"ncp",
"point_sources"
],
"space_time_predictor": {},
"optimised_terms": [],
"optimised_kernel_debugging": [],
"implementation": "generic"
},
"point_sources": 1,
"variables": [
{
"name": "v",
"multiplicity": 2
},
{
"name": "sigma",
"multiplicity": 3
},
{
"name": "lambda",
"multiplicity": 1
},
{
"name": "mu",
"multiplicity": 1
}
],
"plotters": [
{
"type": "vtu::Cartesian::cells::ascii",
"name": "ConservedWriter",
"time": 0.0,
"repeat": 0.05,
"output": "./output/plot",
"variables": 7
},
{
"type": "probe::ascii",
"name": "ProbeWriter",
"time": 0.0,
"repeat": 0.05,
"output": "./output/seismogram-1-",
"variables": 5,
"select": {
"x": 0.2,
"y": 9.99
}
},
{
"type": "probe::ascii",
"name": "ProbeWriter",
"time": 0.0,
"repeat": 0.05,
"output": "./output/seismogram-2-",
"variables": 5,
"select": {
"x": 0.5,
"y": 9.99
}
},
{
"type": "probe::ascii",
"name": "ProbeWriter",
"time": 0.0,
"repeat": 0.05,
"output": "./output/seismogram-3-",
"variables": 5,
"select": {
"x": 0.8,
"y": 9.99
}
}
]
}
]
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "ElasticWaveWriter.h"
ElasticWave::ElasticWaveWriter::ElasticWaveWriter(ElasticWave::MyElasticWaveSolver& solver) {
// @TODO Please insert your code here.
}
ElasticWave::ElasticWaveWriter::~ElasticWaveWriter() {
}
void ElasticWave::ElasticWaveWriter::startPlotting( double time) {
// @TODO Please insert your code here.
}
void ElasticWave::ElasticWaveWriter::finishPlotting() {
// @TODO Please insert your code here.
}
void ElasticWave::ElasticWaveWriter::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* const Q,
double* const 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 was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef POSTPROCESSING_ElasticWaveWriter_CLASS_HEADER_
#define POSTPROCESSING_ElasticWaveWriter_CLASS_HEADER_
#include "exahype/plotters/Plotter.h"
namespace ElasticWave {
class MyElasticWaveSolver;
class ElasticWaveWriter;
}
class ElasticWave::ElasticWaveWriter : public exahype::plotters::Plotter::UserOnTheFlyPostProcessing {
public:
ElasticWaveWriter(ElasticWave::MyElasticWaveSolver& solver);
virtual ~ElasticWaveWriter();
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* const Q,
double* const outputQuantities,
double timeStamp) override;
};
#endif /* POSTPROCESSING_ElasticWaveWriter_CLASS_HEADER_ */
\ No newline at end of file
// This file was generated by the ExaHyPE toolkit.
// It will NOT be regenerated or overwritten.
// Please adapt it to your own needs.
//
// ========================
// www.exahype.eu
// ========================
#include "MyElasticWaveSolver.h"
#include "MyElasticWaveSolver_Variables.h"
#include "kernels/KernelUtils.h"
#include "peano/utils/Loop.h"
tarch::logging::Log ElasticWave::MyElasticWaveSolver::_log( "ElasticWave::MyElasticWaveSolver" );
void ElasticWave::MyElasticWaveSolver::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
initPointSourceLocations(cmdlineargs,constants);
}
void ElasticWave::MyElasticWaveSolver::adjustPointSolution(const double* const x,const double t,const double dt,double* const Q) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
if ( tarch::la::equals( t,0.0 ) ) {
Variables vars(Q);
// Read in values
float lambda = 2.0;
/*std::ifstream fp("x_zero.txt");
if (fp.is_open()) {
while(fp > >string line = ""; getline(infile, line); ){
std::cout << "lambda " << aandb[i] << std::endl;
line >>lambda;
}
}
fp.close();*/
//Initial velocities
Q[0] = 1 * exp(-((x[0] - 5) * (x[0] - 5) + (x[1] - 5) * (x[1] - 5)) / 4);//vx
Q[1] = 1 * exp(-((x[0] - 5) * (x[0] - 5) + (x[1] - 5) * (x[1] - 5)) / 4);//vy
// Initialise stresses to zero
Q[2] = 0.0;//sigma xx
Q[3] = 0.0;//yy
Q[4] = 0.0;//xy
//initialise material parameters
Q[5] = lambda;
Q[6] = 0.5;
}
for(int i=0; i< 7;i++)
assert(std::isfinite(Q[i]));
}
void ElasticWave::MyElasticWaveSolver::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
ReadOnlyVariables varsInside(stateIn);
Variables varsOutside(stateOut);
for (int i = 0; i < (int) sizeof(stateIn); i++){
stateOut[i] = stateIn[i];
}
for (int i = 0; i< (int) sizeof(fluxIn); i++){
fluxOut[i] = fluxIn[i];
}
}
exahype::solvers::Solver::RefinementControl ElasticWave::MyElasticWaveSolver::refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& cellCentre,const tarch::la::Vector<DIMENSIONS,double>& cellSize,double t,const int level) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// Tip: See header file "peano/utils/Loop.h" for dimension-agnostic for loops.
// Example: Loop over all pointwise state variables (plus parameters)
//
// constexpr int sizeOfQ = NumberOfVariables+NumberOfParameters;
// dfor(i,Order+1) {
// const int iLinearised = dLinearised(i,Order+1);
// const double* const Q = luh + iLinearised * sizeOfQ; // pointwise state variables (plus parameters)
// // use Q[0], Q[1], ... Q[sizeOfQ-1]
// }
// @todo Please implement/augment if required
return exahype::solvers::Solver::RefinementControl::Keep;
}
//*****************************************************************************
//******************************** PDE ****************************************
// To use other PDE terms, specify them in the specification file, delete this
// file and its header and rerun the toolkit
//*****************************************************************************
void ElasticWave::MyElasticWaveSolver::eigenvalues(const double* const Q,const int direction,double* const lambda) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
double irho = 1.0; /* density */
double ilambda = Q[5]; /* 1st Lame parameter */
double imu = Q[6]; /* 2nd Lame parameter */
// extract p-wave speed (cp) and s-wave speed (cs)
double cp = std::sqrt((2.0 * imu + ilambda) / irho);
double cs = std::sqrt(imu / irho);
/*
if(normalNonZeroIndex == 0){
lambda[0] = 0.0;
lambda[1] = cp;
lambda[2] = cs;
}else{
lambda[0] = 0.0;
lambda[1] = -cp;
lambda[2] = -cs;
}*/
lambda[0] = cp;
lambda[1] = cs;
lambda[2] = -cp;
lambda[3] = -cs;
lambda[4] = 0.0;
lambda[5] = 0.0;
lambda[6] = 0.0;
}
void ElasticWave::MyElasticWaveSolver::flux(const double* const Q,double** const F) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
Fluxes fluxes(F);
double irho = 1.0; /* density */
double ilambda = Q[5]; /* 1st Lame parameter */
double imu = Q[6]; /* 2nd Lame parameter */
F[0][0] = irho * Q[2];//xx
F[0][1] = irho * Q[4];//xy
F[0][2] = (2.0 * imu + ilambda) * Q[0];
F[0][3] = ilambda * Q[0];
F[0][4] = imu * Q[1];
F[0][5] = 0.0;
F[0][6] = 0.0;
F[1][0] = irho * Q[4];//yx
F[1][1] = irho * Q[3];//yy
F[1][2] = ilambda * Q[1];
F[1][3] = (2.0 * imu + ilambda) * Q[1];
F[1][4] = imu * Q[0];
F[1][5] = 0.0;
F[1][6] = 0.0;
}
void ElasticWave::MyElasticWaveSolver::nonConservativeProduct(const double* const Q,const double* const * const gradQ,double** const BgradQ) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
//Dimensions 2
//Parameters 8
double irho = 1.0; /* density */
double ilambda = Q[5]; /* 1st Lame parameter */
double imu = Q[6]; /* 2nd Lame parameter */
BgradQ[0][0] = 0.0;
BgradQ[0][1] = 0.0;
BgradQ[0][2] = gradQ[0][0];//xx
BgradQ[0][3] = 0.0;//yy
BgradQ[0][4] = gradQ[0][1];//xy
BgradQ[0][5] = 0.0;
BgradQ[0][6] = 0.0;
BgradQ[1][0] = 0.0;
BgradQ[1][1] = 0.0;
BgradQ[1][2] = 0.0;//xx
BgradQ[1][3] = gradQ[1][1];//yy
BgradQ[1][4] = gradQ[1][0];//xy
BgradQ[1][5] = 0.0;
BgradQ[1][6] = 0.0;
}
void ElasticWave::MyElasticWaveSolver::initPointSourceLocations(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants){
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
pointSourceLocation[0][0] = 1.0;
pointSourceLocation[0][1] = 1.0;
}
void ElasticWave::MyElasticWaveSolver::pointSource(const double* const Q,const double* const x,const double t,const double dt, double* const forceVector, int n) {
// Tip: You find documentation for this method in header file "ElasticWave::MyElasticWaveSolver.h".
// Tip: See header file "ElasticWave::AbstractMyElasticWaveSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
forceVector[0] = 1.0;
forceVector[1] = 1.0;
forceVector[2] = 1.0;
forceVector[3] = 1.0;
forceVector[4] = 1.0;
forceVector[5] = 1.0;
forceVector[6] = 1.0;
}
#ifndef __MyElasticWaveSolver_CLASS_HEADER__
#define __MyElasticWaveSolver_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 "AbstractMyElasticWaveSolver.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
namespace ElasticWave{
class MyElasticWaveSolver;
}
class ElasticWave::MyElasticWaveSolver : public ElasticWave::AbstractMyElasticWaveSolver {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
MyElasticWaveSolver(
const double maximumMeshSize,
const int maximumMeshDepth,
const int haloCells,
const int haloBufferCells,
const int limiterBufferCells,
const int regularisedFineGridLevels,
const exahype::solvers::Solver::TimeStepping timeStepping,
const int DMPObservables
);
/**
* 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;
/**
* Adjust the conserved variables and parameters (together: Q) at a given time t at the (quadrature) point x.
*
* @note Please overwrite function adjustSolution(...) if you want to
* adjust the solution degrees of freedom in a cellwise manner.
*
* @param[in] x physical coordinate on the face.
* @param[in] t start of the time interval.
* @param[in] dt width of the time interval.
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
*/
void adjustPointSolution(const double* const x,const double t,const double dt,double* const Q) final override;
/**
* Compute the eigenvalues of the flux tensor per coordinate @p direction.
*
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
* @param[in] direction normal direction of the face / column of the flux vector (range: [0,nDims-1]).
* @param[inout] lambda eigenvalues as C array;
* range: [0,nVar-1], already allocated.
*/
void eigenvalues(const double* const Q,const int direction,double* const lambda) final override;
/**
* Impose boundary conditions at a point on a boundary face
* within the time interval [t,t+dt].
*
* @param[in] x physical coordinate on the face;
* range: [0,nDims-1].
* @param[in] t start of the time interval.
* @param[in] dt 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] direction coordinate direction the face normal is pointing to.
*
* @param[in] QIn the conserved variables (plus parameters) at point x from inside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar+nPar-1], already allocated.
* @param[in] FIn the normal fluxes at point x from inside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar-1], already allocated.
*
* @param[inout] QOut the conserved variables at point x from outside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar+nPar-1], already allocated.
* @param[inout] FOut the normal fluxes at point x from outside of the domain;
* range: [0,nVar-1], already allocated.
*/
void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double * const fluxIn,const double* const stateIn, const double* const gradStateIn,double *fluxOut,double* stateOut) final override;
/**
* Evaluate the refinement criterion within a cell.
*
* @note Instead of a variables array at a single quadrature point we give
* you all NumberOfVariables*(Order+1)^DIMENSIONS solution degrees of freedom.
*
* @note Use this function and ::adjustSolution to set initial conditions.
*
* @param[in] luh all state variables (plus parameters) of a cell;
* range[outer->inner]: [0,(p+1)^nDim-1]x[0,nVar+nPar-1],
* already allocated.
*
* @param[in] cellCentre cellCentre of the cell.
* @param[in] cellSize extent of the cell.
* @param[in] t start of the time interval.
* @param[in] dt width of the time interval.
*
* @return One of exahype::solvers::Solver::RefinementControl::{Erase,Keep,Refine}.
*/
exahype::solvers::Solver::RefinementControl refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& cellCentre,const tarch::la::Vector<DIMENSIONS,double>& cellSize,double t,const int level) override;
//PDE
/**
* Compute the flux tensor.
*
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
*
* @param[inout] F flux at that point;
* range[outer->inner]: [0,nDim-1]x[0,nVar-1],
* already allocated.
*/
void flux(const double* const Q,double** const F) final override;
/* viscousFlux() function not included, as requested in the specification file */
/* algebraicSource() function not included, as requested by the specification file */
/**
* Compute the nonconservative term $B(Q) \nabla Q$.
*
* This function shall return a vector BgradQ which holds the result
* of the full term. To do so, it gets the vector Q and the matrix
* gradQ which holds the derivative of Q in each spatial direction.
* Currently, the gradQ is a continous storage and users can use the
* kernels::icellSize2 class in order to compute the positions inside gradQ.
*
* @TODO: Check if the following is still right:
*
* !!! Warning: BgradQ is a vector of size NumberOfVariables if you
* use the ADER-DG kernels for nonlinear PDEs. If you use
* the kernels for linear PDEs, it is a tensor with dimensions
* Dim x NumberOfVariables.
*
* @param[in] Q vector of state variables (plus material
* parameters); range: [0,nVar+nPar-1], already allocated.
* @param[in] gradQ the gradients of the vector of unknowns, stored
* in a linearized array. (range: [0,dim*(nVar+nPar)-1].
* @param[inout] BgradQ the nonconservative product (extends nVar),
* already allocated.
*/
void nonConservativeProduct(const double* const Q,const double* const * const gradQ,double** const BgradQ) final override;
/**
* Initialise the locations of point sources
* Point source n will be located at pointSourceLocation[n][:]
*
* @param[in] cmdlineargs the command line arguments.
* @param[in] constants access to the constants specified for the solver.
*/
void initPointSourceLocations(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants);
/**
* Compute the contribution of a point source.
*
* This function should return the force vector of a point source n
* located at pointSourceLocation[n][:] as defined in
* initPointSourceLocations()
*
* @param[in] Q vector of state variables (plus material
* parameters), range: [0,nVar+nPar-1], already allocated.
* @param[in] x the position associated with Q
* @param[in] t the time stamp associated with Q
* @param[in] dt the current time step size
* @param[out] forceVector the force vector of the point source (size: ?)
*
* @todo LR: Is x really needed?
* @todo LR: Is dt really needed?
* @todo LR: specify size of forceVector
*/
void pointSource(const double* const Q,const double* const x,const double t,const double dt, double* const forceVector,int n) override;
/* multiplyMaterialParameterMatrix() not included, as requested in the specification file */
};
#endif // __MyElasticWaveSolver_CLASS_HEADER__
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "ProbeWriter.h"
ElasticWave::ProbeWriter::ProbeWriter(ElasticWave::MyElasticWaveSolver& solver) {
// @TODO Please insert your code here.
}
ElasticWave::ProbeWriter::~ProbeWriter() {
}
void ElasticWave::ProbeWriter::startPlotting( double time) {
// @TODO Please insert your code here.
}
void ElasticWave::ProbeWriter::finishPlotting() {
// @TODO Please insert your code here.
}
void ElasticWave::ProbeWriter::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* const Q,
double* const 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 was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef POSTPROCESSING_ProbeWriter_CLASS_HEADER_
#define POSTPROCESSING_ProbeWriter_CLASS_HEADER_
#include "exahype/plotters/Plotter.h"
namespace ElasticWave {
class MyElasticWaveSolver;
class ProbeWriter;
}
class ElasticWave::ProbeWriter : public exahype::plotters::Plotter::UserOnTheFlyPostProcessing {
public:
ProbeWriter(ElasticWave::MyElasticWaveSolver& solver);
virtual ~ProbeWriter();
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* const Q,
double* const outputQuantities,
double timeStamp) override;
};
#endif /* POSTPROCESSING_ProbeWriter_CLASS_HEADER_ */
\ No newline at end of file
export EXAHYPE_CC=mpic++
export COMPILER=GNU
export SHAREDMEM=None
# Make project
export.sh
make
# Compute derivative via finite differencing
./ExaHyPE-ElasticWave ./ElasticWave.exahype2
......@@ -19,8 +19,8 @@
#include "kernels/KernelUtils.h"
#ifdef OPT_KERNELS
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/Kernels.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/converter.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/aderdg/Kernels.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/aderdg/converter.h"
#endif
......@@ -57,6 +57,25 @@ void NavierStokes::NavierStokesSolver_ADERDG::algebraicSource(const tarch::la::V
scenario->source(x, t, ns, Q, S);
}
void NavierStokes::NavierStokesSolver_ADERDG::algebraicSource_vect(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S) {
// naive fallback to scalar
double Q2[NumberOfVariables + NumberOfParameters];
double S2[NumberOfVariables];
for(int i=0; i<VectLength; i++) {
for(int j=0 ; j<NumberOfVariables; j++) {
Q2[j] = Q[j*VectStride+i];
}
for(int j=0; j<NumberOfParameters; j++) {
Q2[j+NumberOfVariables] = P[j*VectStride+i];
}
algebraicSource(x,t,Q2,S2);
for(int j=0 ; j<NumberOfVariables; j++) {
S[j*VectStride+i] = S2[j];
}
}
}
void NavierStokes::NavierStokesSolver_ADERDG::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) {
constexpr auto basisSize = Order + 1;
constexpr auto NumberOfData = NumberOfVariables + NumberOfParameters;
......@@ -385,6 +404,47 @@ void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux(const double *const Q,
ns.evaluateFlux(Q, gradQ, F, true);
}
void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux2(const double *const Q, const double *const P, const double *const *const gradQ, double **F) {
constexpr auto gradQ2Size = NumberOfVariables*DIMENSIONS;
double gradQ2[gradQ2Size];
for(int d = 0; d < DIMENSIONS; d++) {
for(int i=0; i < NumberOfVariables; i++) {
gradQ2[i+NumberOfVariables*d] = gradQ[d][i];
}
}
ns.evaluateFlux(Q, gradQ2, F, true);
}
void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux_vect(const double *const Q, const double *const P, const double *const *const gradQ, double **F) {
// naive fallback to scalar
double Q2[NumberOfVariables + NumberOfParameters];
double gradQ2[NumberOfVariables*DIMENSIONS];
double Fx[NumberOfVariables];
double Fy[NumberOfVariables];
double Fz[NumberOfVariables];
double* F2[3] = {Fx,Fy,Fz};
for(int i=0; i<VectLength; i++) {
for(int j=0 ; j<NumberOfVariables; j++) {
Q2[j] = Q[j*VectStride+i];
}
for(int j=0; j<NumberOfParameters; j++) {
Q2[j+NumberOfVariables] = P[j*VectStride+i];
}
for(int d = 0; d < DIMENSIONS; d++) {
for(int j=0; j < NumberOfVariables; j++) {
gradQ2[j+NumberOfVariables*d] = gradQ[d][j*VectStride+i];
}
}
viscousFlux(Q2,gradQ2,F2);
for(int d=0 ; d<DIMENSIONS; d++) {
for(int j=0; j < NumberOfVariables; j++) {
F[d][j*VectStride+i] = F2[d][j];
}
}
}
}
void NavierStokes::NavierStokesSolver_ADERDG::boundaryConditions( double* const fluxIn, const double* const stateIn, const double* const gradStateIn, const double* const luh, const tarch::la::Vector<DIMENSIONS, double>& cellCentre, const tarch::la::Vector<DIMENSIONS,double>& cellSize, const double t,const double dt, const int direction, const int orientation) {
#ifdef OPT_KERNELS
......
......@@ -156,9 +156,13 @@ class NavierStokes::NavierStokesSolver_ADERDG : public NavierStokes::AbstractNav
*/
//void flux(const double* const Q,double** F) final override;
void viscousFlux(const double* const Q, const double* const gradQ, double** F) final override;
void viscousFlux2(const double* const Q, const double* const P, const double* const * const gradQ, double** F);
void viscousFlux_vect(const double* const Q, const double* const P, const double* const * const gradQ, double** F);
// TODO(Lukas) Add to toolkit!
void algebraicSource(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, double *S) override;
inline void algebraicSource2(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S) {algebraicSource(x,t,Q,S);}
void algebraicSource_vect(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S);
/* nonConservativeProduct() function is not included, as requested in the specification file */
......
#ifndef __CurvilinearTransformation_CLASS_HEADER__
#define __CurvilinearTransformation_CLASS_HEADER__
#include <array>
#include "tarch/la/Vector.h"
#include "kernels/aderdg/generic/Kernels.h"
struct Boundary_single_coordinate{
std::vector<double> left;
std::vector<double> right;
std::vector<double> bottom;
std::vector<double> top;
std::vector<double> front;
std::vector<double> back;
Boundary_single_coordinate(int nx, int ny, int nz){
left .resize(ny*nz);
right .resize(ny*nz);
bottom.resize(nx*nz);
top .resize(nx*nz);
front .resize(nx*ny);
back .resize(nx*ny);
}
~Boundary_single_coordinate(){
}
friend std::ostream& operator<< (std::ostream& stream, const Boundary_single_coordinate& boundary){
stream << std::endl;
stream << "Left" << std::endl;
for(int i = 0 ; i<boundary.left.size(); i++){
stream << boundary.left[i] << std::endl;
}
stream << "Right" << std::endl;
for(int i = 0 ; i<boundary.right.size(); i++){
stream << boundary.right[i] << std::endl;
}
stream << "Bottom" << std::endl;
for(int i = 0 ; i<boundary.bottom.size(); i++){
stream << boundary.bottom[i] << std::endl;
}
stream << "Top" << std::endl;
for(int i = 0 ; i<boundary.top.size(); i++){
stream << boundary.top[i] << std::endl;
}
stream << "Front" << std::endl;
for(int i = 0 ; i<boundary.front.size(); i++){
stream << boundary.front[i] << std::endl;
}
stream << "Back" << std::endl;
for(int i = 0 ; i<boundary.back.size()-1; i++){
stream << boundary.back[i] << std::endl;
}
return stream << boundary.back[boundary.back.size()-1] << std::endl;
}
};
class CurvilinearTransformation{
public:
//Default Constructor
CurvilinearTransformation():_num_nodes(0), _dx(0),
_fault_position(0),
_right_vertex{0.},
_left_vertex{0.}{};
CurvilinearTransformation(const int num_nodes, const int mesh_level ,
const double fault_position,
double* domain_offset,
double* domain_size);
//Destructor
~CurvilinearTransformation(){
delete[] _boundary_x;
delete[] _boundary_y;
delete[] _boundary_z;
};
void genCoordinates(const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
double* gl_vals_x,double* gl_vals_y,double* gl_vals_z,
double* jacobian,
double* q_x,double* q_y,double* q_z,
double* r_x,double* r_y,double* r_z,
double* s_x,double* s_y,double* s_z);
void metricDerivativesAndJacobian3D(const double* const curvilinear_x,
const double* const curvilinear_y,
const double* const curvilinear_z,
const double* const dx,
double* gl_vals_x, double* gl_vals_y, double* gl_vals_z,
double* q_x, double* q_y, double* q_z,
double* r_x, double* r_y, double* r_z,
double* s_x, double* s_y, double* s_z,
double* jacobian);
private:
int _num_nodes;
double _dx;
double _block_width_x[2];
double _fault_position;
std::array<double,DIMENSIONS> _left_vertex;
std::array<double,DIMENSIONS> _right_vertex;
std::array<double,DIMENSIONS> _rect_width;
//number of elements over whole domain
int _ne_x_block[2];
int _ne_x;
int _ne_y;
int _ne_z;
//number of nodes over whole domain
int _nx_block[2];
int _nx;
int _ny;
int _nz;
//boundaries
Boundary_single_coordinate* _boundary_x[2];
Boundary_single_coordinate* _boundary_y[2];
Boundary_single_coordinate* _boundary_z[2];
double* lagrange_basis_at_nodes;
double* denominator_lagrange;
double* unif_mesh;
double fault(double y, double z);
double topography(double x, double z, double depth);
/* #if defined(USE_ASAGI) */
/* double topography_fromASAGI(double x, double z, double* topography, easi::ArraysAdapter& adapter, easi::Component* model); */
/* #endif */
void getInterpolatedFace_fromBottomAndTop( int n_block, int n_left_right, int n_top_bottom, int face,
double* top_bnd_x, double* top_bnd_y, double* top_bnd_z,
double* bottom_bnd_x, double* bottom_bnd_y, double* bottom_bnd_z,
double* face_x, double* face_y , double* face_z);
double interpolate2D(double x, double y,
int number_nodes_x, int number_nodes_y,
const double* const nodes_x, const double* const nodes_y,
const double* const values);
double interpolate2D_uniform(double x, double y, double* values);
double interpolate3D_uniform(int x, int y ,int z, const double* const values);
double lagrangeBasis(double x, int i, int n, const double* const nodes);
double lagrangeBasis_uniform(double x,int i);
void getValuesAtQuadNodes2D(double* dest_mesh, double* results);
void getValuesAtQuadNodes3D(const double* const dest_mesh, double* results);
void computeDerivatives_x(int i, int j, double* values , int num_nodes, double& der_x, double dx);
void computeDerivatives_y(int i, int j, double* values , int num_nodes, double& der_y, double dy);
void computeDerivatives_x_3D(int i, int j, int k, double* values , int num_nodes, double& der_x, double dx);
void computeDerivatives_y_3D(int i, int j, int k, double* values , int num_nodes, double& der_y, double dy);
void computeDerivatives_z_3D(int i, int j, int k ,double* values , int num_nodes, double& der_z, double dz);
void metricDerivativesAndJacobian(int num_nodes, double* curvilinear_x, double* curvilinear_y,double* gl_vals_x,double* gl_vals_y,double* q_x,double* q_y,double* r_x,double* r_y,double* jacobian, double dx, double dy);
double interpolate_fault_surface(const double* const fault_x,
const double* const fault_y,
const double* const fault_z,
const double y, const double z);
double interpol2d_dG(int my, int mz, int nmy, int npy, int nmz, int npz, double y, double z, double* yj, double* zk, double* f);
double lagrange(int m, int p, int i, double x, double *xi);
void getBoundaryCurves3D_cutOffTopography_withFault(int n,
Boundary_single_coordinate* boundary_x,
Boundary_single_coordinate* boundary_y,
Boundary_single_coordinate* boundary_z);
void transFiniteInterpolation(int mx, int my, int j_m, int j_p, int i_m, int i_p, int num_points, double* left_bnd_x, double* left_bnd_y, double* right_bnd_x, double* right_bnd_y, double* bottom_bnd_x, double* bottom_bnd_y, double* top_bnd_x, double* top_bnd_y, double* curvilinear_x , double* curvilinear_y );
void transFiniteInterpolation_singleCoordinate(int mx, int my,
double* left_bnd, double* right_bnd,
double* bottom_bnd, double* top_bnd,
double* curvilinear);
void transFiniteInterpolation_singleCoordinate(int mx, int my,
std::vector<double> left_bnd,
std::vector<double> right_bnd,
std::vector<double> bottom_bnd,
std::vector<double> top_bnd,
double* curvilinear){
transFiniteInterpolation_singleCoordinate(mx, my,
left_bnd.data(), right_bnd.data(),
bottom_bnd.data(), top_bnd.data(),
curvilinear);
};
void transFiniteInterpolation3D(int n,
int k_m, int k_p ,
int j_m, int j_p ,
int i_m, int i_p ,
Boundary_single_coordinate* boundary,
double* curvilinear);
int getBlock(const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx);
};
#endif
/**
* 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
**/
/**
Elastic Wave
A simple project. (well, it was simple, in the beginning).
*/
exahype-project Elastic
peano-kernel-path const = ./Peano
exahype-path const = ./ExaHyPE
output-directory const = ./ApplicationExamples/Linear/ElasticWave2D
architecture const = noarch
log-file = mylogfile.log
computational-domain
dimension const = 2
width = 10.0, 10.0
offset = 0.0, 0.0
end-time = 10.0
end computational-domain
/*shared-memory
identifier = dummy
configure = {}
cores = 8
properties-file = sharedmemory.properties
end shared-memory*/
/*distributed-memory
identifier = static_load_balancing
configure = {greedy-naive,FCFS}
// configure = {hotspot,fair,ranks-per-node:4}
buffer-size = 64
timeout = 60
end distributed-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 ADER-DG MyElasticWaveSolver
variables const = velocity:2,stress:3
parameters const = density:1,pwavespeed:1,swavespeed:1
order const = 3
/* 27 points: 0.05, 9 points: 0.15 */
maximum-mesh-size = 0.5
maximum-mesh-depth = 0
time-stepping = global
type const = linear, Legendre
terms const = flux,ncp,materialparameters,pointsources:2
optimisation const = generic,patchwiseadjust
language const = C
//plot Peano::Legendre::cells::hdf5 ConservedQuantitiesWriter
//plot Peano::Legendre::cells::ascii ConservedQuantitiesWriter
//plot vtk::Legendre::cells::binary ConservedQuantitiesWriter
plot vtk::Legendre::cells::ascii ConservedQuantitiesWriter
variables const = 3
time = 0.0
repeat = 1e-2
output = ./conserved
end plot
/* this is the fake plotter used to compute global integrals */
/* it has no output fields. */
plot vtk::Cartesian::vertices::ascii ComputeGlobalIntegrals
variables const = 0
time = 0.0
repeat = 0.05
output = ./output/these-files-should-not-be-there
end plot
plot vtk::Cartesian::vertices::ascii PrimitivesWriter
variables const = 3
time = 0.0
repeat = 0.05
output = ./primitive
end plot
plot vtk::Cartesian::vertices::ascii ExactPrimitivesWriter
variables const = 3
time = 0.0
repeat = 0.05
output = ./exact-primitive
end plot
/* Do not need the time series for a point in the moment
plot probe::ascii ProbeWriter
variables const = 5
time = 0.0
repeat = 0.05
output = ./seismogram
select = x:0.2,y:0.2
end plot
*/
end solver
end exahype-project
#ifndef __MyElasticWaveSolver_CLASS_HEADER__
#define __MyElasticWaveSolver_CLASS_HEADER__
// 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
// ========================
#include <ostream>
#include "AbstractMyElasticWaveSolver.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
#include "tarch/la/Vector.h"
namespace Elastic{
class MyElasticWaveSolver;
}
class Elastic::MyElasticWaveSolver : public Elastic::AbstractMyElasticWaveSolver {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
MyElasticWaveSolver(const double maximumMeshSize,const int maximumMeshDepth,const int haloCells,const int regularisedFineGridLevels,const exahype::solvers::Solver::TimeStepping timeStepping,const int DMPObservables);
/**
* Initialise the solver.
*
* \param[in] cmdlineargs the command line arguments.
*/
void init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) final override;
/**
* Patchwise adjust
* @TODO LR : Document
*/
void adjustSolution(double* const luh,const tarch::la::Vector<DIMENSIONS,double>& center,const tarch::la::Vector<DIMENSIONS,double>& dx,double t,double dt) final 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) final 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[in] FIn the normal fluxes 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).
* \param[inout] FOut the normal fluxes 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) final override;
/**
* Evaluate the refinement criterion within a cell.
*
* \note Instead of a variables array at a single quadrature point we give
* you all NumberOfVariables*(Order+1)^DIMENSIONS solution degrees of freedom.
*
* \note Use this function and ::adjustSolution to set initial conditions.
*
* \param[in] centre The centre of the cell.
* \param[in] dx The extent of the cell.
* \param[in] t the start of the time interval.
* \param[in] dt the width of the time interval.
* \return One of exahype::solvers::Solver::RefinementControl::{Erase,Keep,Refine}.
*/
exahype::solvers::Solver::RefinementControl refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& centre,const tarch::la::Vector<DIMENSIONS,double>& dx,double t,const int level) override;
//PDE
/**
* Compute the flux tensor.
*