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

Commit 3d2bafda authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard

Merge branch 'master' into jm/splitPicard

parents 9ffeff65 15bf23fe
// 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.
//
//
// ========================