Commit b2f4e388 authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard
Browse files

ApplicationExamples - add bubble vect test

parent fceada1f
{
"project_name": "NavierStokes",
"paths": {
"peano_kernel_path": "./Peano/",
"exahype_path": "./ExaHyPE",
"output_directory": "./ApplicationExamples/TestCases/nonlinear_vectorized_navier_bubble_3d"
},
"computational_domain": {
"dimension": 3,
"time_steps": 200,
"offset": [
0.0,
0.0,
0.0
],
"width": [
1000.0,
1000.0,
1000.0
]
},
"architecture": "hsw",
"shared_memory": {
"cores": 12,
"background_job_consumers": 12,
"properties_file": "sharedmemory.properties",
"autotuning_strategy": "dummy"
},
"_distributed_memory": {
"ranks_per_node": 4,
"load_balancing_strategy": "hotspot",
"load_balancing_type": "static",
"node_pool_strategy": "fair",
"timeout": 60000,
"buffer_size": 64
},
"optimisation": {
"fuse_algorithmic_steps": "all",
"fuse_algorithmic_steps_diffusion_factor": 0.99,
"fuse_algorithmic_steps_rerun_factor": 0.99,
"spawn_predictor_as_background_thread": false,
"spawn_amr_background_threads": true,
"spawn_update_as_background_thread": true,
"time_step_batch_factor": 1.0,
"disable_metadata_exchange_in_batched_time_steps": true,
"disable_vertex_exchange_in_time_steps": true,
"double_compression": 0.0,
"spawn_double_compression_as_background_thread": false
},
"solvers": [
{
"type": "ADER-DG",
"name": "NavierStokesSolver_ADERDG",
"order": 6,
"maximum_mesh_size": 300,
"maximum_mesh_depth": 0,
"time_stepping": "globalfixed",
"halo_cells": 0,
"aderdg_kernel": {
"language": "C",
"nonlinear": true,
"terms": [
"flux",
"viscous_flux",
"source"
],
"space_time_predictor": {"predictor_recompute": true, "vectorise_terms":true, "AoSoA2_layout":true},
"optimised_terms": [],
"optimised_kernel_debugging": ["flops"],
"implementation": "optimised"
},
"point_sources": 0,
"variables": [
{
"name": "rho",
"multiplicity": 1
},
{
"name": "j",
"multiplicity": 3
},
{
"name": "E",
"multiplicity": 1
}
],
"material_parameters": [
{
"name": "height",
"multiplicity": 1
},
{
"name": "backgroundstate",
"multiplicity": 2
}
],
"parameters": {
"use-gravity": true,
"use-background-state": true,
"viscosity": 0.01,
"scenario": "two-bubbles"
},
"_plotters": [
{
"type": "vtu::Cartesian::vertices::ascii",
"name": "Plotter",
"time": 10.0,
"repeat": 5.0,
"output": "plots/solution_cartesian",
"variables": 6
},
{
"type": "vtu::Legendre::vertices::ascii",
"name": "Plotter",
"time": 10.0,
"repeat": 1.0,
"output": "plots/solution",
"variables": 6
}
]
}
]
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "ErrorWriter.h"
#include "kernels/GaussLegendreBasis.h"
#include "kernels/KernelUtils.h"
#include "peano/utils/Loop.h"
#include "tarch/la/VectorOperations.h"
#include <algorithm>
#include <fstream>
#include <iomanip>
#include <sstream>
#include "PDE.h"
NavierStokes::ErrorWriter::ErrorWriter(NavierStokesSolver_ADERDG& solver)
: solver(&solver),
exahype::plotters::ADERDG2UserDefined::ADERDG2UserDefined(),
hmin(std::numeric_limits<double>::max()) {}
void NavierStokes::ErrorWriter::plotPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* u,
double timeStamp) {
if (solver->scenarioName == "two-bubbles") {
// Only measure time!
return;
}
constexpr int numberOfVariables = NavierStokesSolver_ADERDG::NumberOfVariables;
constexpr int numberOfParameters = NavierStokesSolver_ADERDG::NumberOfParameters;
constexpr int numberOfData = numberOfVariables + numberOfParameters;
constexpr int order = NavierStokesSolver_ADERDG::Order;
constexpr int basisSize = order + 1;
constexpr int quadOrder = 9;
constexpr int numberOfQuadPoints = quadOrder + 1;//basisSize;
constexpr int gradSize = DIMENSIONS * numberOfVariables;
const auto& quadratureNodes = kernels::legendre::nodes[quadOrder];
const auto& quadratureWeights = kernels::legendre::weights[quadOrder];
assertionEqualsMsg(DIMENSIONS, 2, "ErrorWriter is only supported for 2D!");
double x[2] = {0.0, 0.0};
kernels::idx4 idx(basisSize, basisSize, basisSize, numberOfData);
dfor(i, numberOfQuadPoints) {
long double w_dV = 1.0;
for (int d = 0; d < DIMENSIONS; d++) {
x[d] = offsetOfPatch[d] +
sizeOfPatch[d] * quadratureNodes[i(d)];
w_dV *= sizeOfPatch[d] * quadratureWeights[i(d)];
hmin =
std::min(hmin, sizeOfPatch[d]); // TODO(Lukas) Is this what we want?
}
// Compute analytical solution
auto uAna = std::array<double, numberOfVariables>{};
auto uAnaGrad = std::array<double, gradSize>{};
auto vars = Variables(uAna.data());
solver->scenario->analyticalSolution(x, timeStamp, solver->ns, vars, uAnaGrad.data());
for (int v = 0; v < numberOfVariables; v++) {
const auto curAna = uAna[v];
// TODO(Lukas) Correct interpolation?
const auto curNum = kernels::legendre::interpolate(offsetOfPatch.data(),
sizeOfPatch.data(),
x,
numberOfData,
v,
order,
u);
const long double uDiff = std::abs(
static_cast<long double>(curNum) -
static_cast<long double>(curAna));
errorL2[v] += uDiff * uDiff * w_dV;
errorL1[v] += uDiff * w_dV;
errorLInf[v] = std::max(errorLInf[v], uDiff);
normL1Ana[v] += std::abs(uAna[v]) * w_dV;
normL2Ana[v] += uAna[v] * uAna[v] * w_dV;
normLInfAna[v] = std::max(normLInfAna[v],
static_cast<long double>(std::abs(uAna[v])));
}
}
}
void NavierStokes::ErrorWriter::startPlotting(double time) {
timeStamp = time;
std::fill(errorL1.begin(), errorL1.end(), 0.0);
std::fill(errorL2.begin(), errorL2.end(), 0.0);
std::fill(errorLInf.begin(), errorLInf.end(), 0.0);
std::fill(normL1Ana.begin(), normL1Ana.end(), 0.0);
std::fill(normL2Ana.begin(), normL2Ana.end(), 0.0);
std::fill(normLInfAna.begin(), normLInfAna.end(), 0.0);
}
void NavierStokes::ErrorWriter::finishPlotting() {
constexpr int numberOfVariables =
AbstractNavierStokesSolver_ADERDG::NumberOfVariables;
// Compute elapsed simulation time.
const auto endTime = std::chrono::high_resolution_clock::now();
const auto timeDiff = std::chrono::duration_cast<std::chrono::milliseconds>(
endTime - startTime).count();
// Check whether our file exists already.
const bool isExisting = [&]() -> bool {
auto file = std::ifstream(filename);
return file.good();
}();
auto file = std::ofstream(filename, std::ios::app);
if (!file.is_open()) {
throw std::runtime_error("ErrorWriter: Error while opening file " +
filename + "!");
}
assert(file.is_open());
// TODO(Lukas) Is this enough precision?
const auto precision = std::numeric_limits<long double>::max_digits10 + 5;
// Write csv header in case of new file.
if (!isExisting) {
file << std::setprecision(precision) << std::fixed;
file << "nbackgroundJobs,norm,order,hmin,";
for (int i = 0; i < numberOfVariables; ++i) {
file << "var" << i << ",";
}
file << "time" << "," << "runTime" << std::endl;
}
auto plotRow = [&](const std::string& normType, const Array_t arr) {
file << std::setprecision(precision) << std::fixed
<< solver->MaxNumberOfRunningBackgroundJobConsumerTasksDuringTraversal << ","
<< normType << ","
<< AbstractNavierStokesSolver_ADERDG::Order << "," << hmin << ",";
for (int i = 0; i < numberOfVariables; ++i) {
file << std::setprecision(precision) << std::fixed << arr[i] << ",";
}
file << timeStamp << "," << timeDiff << std::endl;
};
// Note that we write the l2-norms without applying the square root.
// The reason is expected precision loss when using MPI.
// (need to aggregate results of multiple nodes!)
plotRow("l1", errorL1);
plotRow("l2", errorL2);
plotRow("lInf", errorLInf);
plotRow("l1Norm", normL1Ana);
plotRow("l2Norm", normL2Ana);
plotRow("lInfNorm", normLInfAna);
// Verify that we actually wrote something.
file << std::flush; // Buffering might hide error
if (!file.good()) {
throw std::runtime_error("ErrorWriter: Error while writing to file " +
filename + "!");
}
// Now pretty print for console:
auto prettyPrintRow = [&](const std::string& normType, const Array_t arr) {
std::cout << "[" << timeDiff << "ms]: " << normType << "-Error for order "
<< AbstractNavierStokesSolver_ADERDG::Order << " and hmin of " << hmin
<< " at time " << timeStamp << ":" << std::endl;
for (int i = 0; i < numberOfVariables; ++i) {
std::cout << std::setprecision(5) << std::scientific << arr[i] << "\t";
}
std::cout << std::endl;
};
// Compute square roots for pretty printing.
for (int v = 0; v < numberOfVariables; v++) {
errorL2[v] = sqrt(errorL2[v]);
normL2Ana[v] = sqrt(normL2Ana[v]);
}
std::cout << "\n\n";
prettyPrintRow("l1", errorL1);
prettyPrintRow("l2", errorL2);
prettyPrintRow("lInf", errorLInf);
prettyPrintRow("l1Norm", normL1Ana);
prettyPrintRow("l2Norm", normL2Ana);
prettyPrintRow("lInfNorm", normLInfAna);
std::cout << "\n\n";
}
void NavierStokes::ErrorWriter::init(
const std::string& filename, int orderPlusOne, int solverUnknowns,
int writtenUnknowns, exahype::parser::ParserView plotterParameters) {
ADERDG2UserDefined::init(filename, orderPlusOne, solverUnknowns,
writtenUnknowns, plotterParameters);
const auto& node = tarch::parallel::Node::getInstance();
const auto rank = node.getRank();
const auto no_nodes = node.getNumberOfNodes();
isMpi = no_nodes > 1;
auto fs = std::stringstream();
fs << filename;
if (isMpi) {
fs << "_rank_" << rank;
}
fs << ".csv";
this->filename = fs.str();
startTime = std::chrono::high_resolution_clock::now();
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef ErrorWriter_CLASS_HEADER_
#define ErrorWriter_CLASS_HEADER_
#include "NavierStokesSolver_ADERDG.h"
#include "exahype/plotters/ADERDG2UserDefined.h"
#include <array>
#include <chrono>
namespace NavierStokes {
class ErrorWriter;
}
class NavierStokes::ErrorWriter : public exahype::plotters::ADERDG2UserDefined {
private:
double timeStamp;
std::string filename;
bool isMpi;
NavierStokesSolver_ADERDG* solver;
using Array_t = std::array<long double, NavierStokesSolver_ADERDG::NumberOfVariables>;
double hmin;
Array_t errorL1, errorL2, errorLInf;
Array_t normL1Ana, normL2Ana, normLInfAna;
using Time_t = std::chrono::time_point<std::chrono::high_resolution_clock>;
Time_t startTime;
public:
void init(const std::string& filename, int orderPlusOne, int solverUnknowns,
int writtenUnknowns,
exahype::parser::ParserView plotterParameters) override;
/**
* Constructor.
*
* \note ExaHyPE does not increment file counters for
* you if you use user defined plotting. You have
* to declare and manage such member variables yourself.
*/
ErrorWriter(NavierStokesSolver_ADERDG& solver);
/**
* This method is invoked every time a cell
* is touched by the plotting device.
*
* \note Use the protected variables _order, _variables to
* determine the size of u.
* The array u has the size _variables * (_order+1)^DIMENSIONS.
*
* \param[in] offsetOfPatch the offset of the cell/patch.
* \param[in] sizeOfPatch the offset of the cell/patch.
* \param[in] u the degrees of freedom "living" inside of the patch.
*/
void plotPatch(const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
double* u, double timeStamp) override;
/**
* This method is called at the beginning of the plotting.
* You can use it to reset member variables, e.g., those
* used for calculations, or to increment file counters.
*
* \param[in] time a characteristic solver time stamp.
* Usually the global minimum.
*/
void startPlotting(double time) override;
/**
* This method is called at the end of the plotting.
* You can use it to reset member variables, finalise calculations (compute
* square roots etc.), or to increment file counters
*/
void finishPlotting() override;
};
#endif /* ErrorWriter_CLASS_HEADER_ */
#ifndef __NavierStokesSolverDG_CLASS_HEADER__
#define __NavierStokesSolverDG_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 "AbstractNavierStokesSolver_ADERDG.h"
#include <string>
#include "Scenarios/Scenario.h"
#include "kernels/GaussLegendreBasis.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include <memory>
#include "tarch/logging/Log.h"
#include "PDE.h"
//#include "AMR/AMRSettings.h"
namespace NavierStokes{
class NavierStokesSolver_ADERDG;
}
class NavierStokes::NavierStokesSolver_ADERDG : public NavierStokes::AbstractNavierStokesSolver_ADERDG {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
static double** weights;
std::unique_ptr<Scenario> scenario;
std::string scenarioName;
PDE ns;
//AMRSettings amrSettings;
NavierStokesSolver_ADERDG(
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 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[inout] Q the conserved variables (and parameters) associated with a quadrature point
* as C array (already allocated).
*/
void adjustPointSolution(const double* const x,const double t,const double dt,double* Q) 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* lambda) final override;
void viscousEigenvalues(const double* const Q,const int d,double* 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;
void 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);
void mapDiscreteMaximumPrincipleObservables(double* observables,const double* const Q) const override;
bool isPhysicallyAdmissible(
const double* const solution,
const double* const observablesMin,const double* const observablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
const double t) const override;
/**