04.06., 9:00 - 12:00: GitLab will be migrated to a new server environment and upgraded to Enterprise Edition Ultimate. The estimated downtime will be 2-3 hours. Please see https://doku.lrz.de/display/PUBLIC/GitLab+Ultimate+Migration for more details about changes related to the migration.

Commit 163e9dbf authored by Anne Reinarz's avatar Anne Reinarz

add reader

parent 9344d734
#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.
*
* \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) final override;
/**
* 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::idx2 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 the vector of unknowns at the given position
* \param[in] gradQ the gradients of the vector of unknowns,
* stored in a linearized array.
* \param[inout] The vector BgradQ (extends nVar), already allocated.
*
**/
void nonConservativeProduct(const double* const Q,const double* const gradQ,double* const BgradQ) final override;
/* void coefficientMatrix(const double* const Q,const int d,double* const Bn); */
void pointSource(const double* const x,const double t,const double dt, double* const forceVector, double* const x0, int n) final override;
void riemannSolver(double* const FL,double* const FR,const double* const QL,const double* const QR,const double t,const double dt, const tarch::la::Vector<DIMENSIONS, double>& lengthScale,const int direction, bool isBoundaryFace, int faceIndex) override;
/**
* @TODO LR : document
*/
void multiplyMaterialParameterMatrix(const double* const Q, double* const rhs) final override;
/* implement user defined numerical boundary procedures*/
void riemannSolver_BC0(double v, double sigma, double z, double r, double& v_hat, double& sigma_hat);
void riemannSolver_BCn(double v, double sigma, double z, double r, double& v_hat, double& sigma_hat);
/* implement user defined numerical fluxes*/
void riemannSolver_Nodal(double v_p,double v_m, double sigma_p, double sigma_m, double z_p , double z_m, double& v_hat_p , double& v_hat_m, double& sigma_hat_p, double& sigma_hat_m);
void localBasis(double* const n, double* const m);
void riemannSolver_boundary(int faceIndex,double r, double vn , double vm, double Tn , double Tm, double zp, double zs, double& vn_hat , double& vm_hat , double& Tn_hat , double& Tm_hat);
void generate_fluctuations_right(double z, double T,double T_hat,double v, double v_hat, double& F);
void generate_fluctuations_left(double z, double T,double T_hat,double v, double v_hat, double& F);
void rotate_into_physical_basis(double* const n,double* const m, double Fn,double Fm, double& Fx, double& Fy);
void rotate_into_orthogonal_basis(double* const n,double* const m, double Tx,double Ty, double& Tn, double& Tm);
void extract_tractions_and_particle_velocity(double* const n, const double* const Q, double& Tx,double& Ty,double& vx,double& vy);
void get_normals(int normalNonZeroIndex,double& norm, double* const n,const double* const Q);
};
#endif // __MyElasticWaveSolver_CLASS_HEADER__
/**
* 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_Topo
architecture const = noarch
log-file = mylogfile.log
computational-domain
dimension const = 2
width = 1.0, 1.0
offset = 0.0, 0.0
end-time = 4.0
end computational-domain
shared-memory
identifier = dummy
configure = {}
cores = 4
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 = v:2,sigma:3
parameters const = rho:1,cp:1,cs:1,jacobian:1,metric_derivative:4,curve_grid:2
order const = 4
/* 27 points: 0.05, 9 points: 0.15 */
maximum-mesh-size = 0.05
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::vertices::ascii ConservedQuantitiesWriter
variables const = 15
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
#include "kernels/KernelUtils.h"
#include "kernels/GaussLegendreQuadrature.h"
//namespace Linear{
// class CurvilinearTransformation;
//}
//class Linear::CurvilinearTransformation {