C++ manager class for PdeFiniteDifferenceKernels API. The low level calls are managed in the namespace pde::detail DeviceManager, whereas the high level infrastructure is delegated to the particular solver type.
Only linear hyperbolic and parabolic PDEs are supported (up to 3D). The exposed implementation is through:
- AdvectionDiffusionSolver1D, AdvectionDiffusionSolver2D
- WaveEquationSolver1D, WaveEquationSolver2D
These solvers are implemented with the Curiously Recurring Template Pattern (CRTP), useful for delegating the members data type at compile time.
For convenience's sake the following typedefs have been defined:
- Advection-Diffusion:
typedef AdvectionDiffusionSolver1D<MemorySpace::Device, MathDomain::Float> GpuSingleAdvectionDiffusionSolver1D;
typedef GpuSingleAdvectionDiffusionSolver1D GpuFloatAdvectionDiffusionSolver1D;
typedef AdvectionDiffusionSolver1D<MemorySpace::Device, MathDomain::Double> GpuDoubleAdvectionDiffusionSolver1D;
typedef AdvectionDiffusionSolver1D<MemorySpace::Host, MathDomain::Float> CpuSingleAdvectionDiffusionSolver1D;
typedef CpuSingleAdvectionDiffusionSolver1D CpuFloatAdvectionDiffusionSolver1D;
typedef AdvectionDiffusionSolver1D<MemorySpace::Host, MathDomain::Double> CpuDoubleAdvectionDiffusionSolver1D;
typedef GpuSingleAdvectionDiffusionSolver1D ad1D;
typedef GpuDoubleAdvectionDiffusionSolver1D dad1D;
typedef AdvectionDiffusionSolver2D<MemorySpace::Device, MathDomain::Float> GpuSingleAdvectionDiffusionSolver2D;
typedef GpuSingleAdvectionDiffusionSolver2D GpuFloatAdvectionDiffusionSolver2D;
typedef AdvectionDiffusionSolver2D<MemorySpace::Device, MathDomain::Double> GpuDoubleAdvectionDiffusionSolver2D;
typedef AdvectionDiffusionSolver2D<MemorySpace::Host, MathDomain::Float> CpuSingleAdvectionDiffusionSolver2D;
typedef CpuSingleAdvectionDiffusionSolver2D CpuFloatAdvectionDiffusionSolver2D;
typedef AdvectionDiffusionSolver2D<MemorySpace::Host, MathDomain::Double> CpuDoubleAdvectionDiffusionSolver2D;
typedef GpuSingleAdvectionDiffusionSolver2D ad2D;
typedef GpuDoubleAdvectionDiffusionSolver2D dad2D;- Wave Equation:
typedef WaveEquationSolver1D<MemorySpace::Device, MathDomain::Float> GpuSingleWaveEquationSolver1D;
typedef GpuSingleWaveEquationSolver1D GpuFloatWaveEquationSolver1D;
typedef WaveEquationSolver1D<MemorySpace::Device, MathDomain::Double> GpuDoubleWaveEquationSolver1D;
typedef WaveEquationSolver1D<MemorySpace::Host, MathDomain::Float> CpuSingleWaveEquationSolver1D;
typedef CpuSingleWaveEquationSolver1D CpuFloatWaveEquationSolver1D;
typedef WaveEquationSolver1D<MemorySpace::Host, MathDomain::Double> CpuDoubleSolver1D;
typedef GpuSingleWaveEquationSolver1D wave1D;
typedef GpuDoubleWaveEquationSolver1D dwave1D;
typedef WaveEquationSolver2D<MemorySpace::Device, MathDomain::Float> GpuSingleWaveEquationSolver2D;
typedef GpuSingleWaveEquationSolver2D GpuFloatWaveEquationSolver2D;
typedef WaveEquationSolver2D<MemorySpace::Device, MathDomain::Double> GpuDoubleWaveEquationSolver2D;
typedef WaveEquationSolver2D<MemorySpace::Host, MathDomain::Float> CpuSingleWaveEquationSolver2D;
typedef CpuSingleWaveEquationSolver2D CpuFloatWaveEquationSolver2D;
typedef WaveEquationSolver2D<MemorySpace::Host, MathDomain::Double> CpuDoubleSolver2D;
typedef GpuSingleWaveEquationSolver2D wave2D;
typedef GpuDoubleWaveEquationSolver2D dwave2D; cl::vec grid = cl::LinSpace(0.0f, 1.0f, 128);
auto _grid = grid.Get();
std::vector<float> _initialCondition(grid.size());
for (unsigned i = 0; i < _initialCondition.size(); ++i)
_initialCondition[i] = sin(_grid[i]);
cl::vec initialCondition(_initialCondition);
unsigned steps = 100;
double dt = 1e-4;
float velocity = .05f;
float diffusion = 0.1f;
BoundaryCondition leftBoundaryCondition(BoundaryConditionType::Neumann, 0.0);
BoundaryCondition rightBoundaryCondition(BoundaryConditionType::Neumann, 0.0);
BoundaryCondition1D boundaryConditions(leftBoundaryCondition, rightBoundaryCondition);
pde::GpuSinglePdeInputData1D data(initialCondition, grid, velocity, diffusion, dt, solverType, SpaceDiscretizerType::Centered, boundaryConditions);
pde::ad1D solver(data);
solver.Advance(steps);
const auto solution = solver.solution->columns[0]->Get(); cl::dvec grid = cl::LinSpace<MemorySpace::Device, MathDomain::Double>(0.0, 1.0, 128);
auto _grid = grid.Get();
std::vector<double> _initialCondition(grid.size());
for (unsigned i = 0; i < _initialCondition.size(); ++i)
_initialCondition[i] = sin(_grid[i]);
cl::dvec initialCondition(_initialCondition);
unsigned steps = 100;
double dt = 1e-4;
double velocity = .05;
BoundaryCondition leftBoundaryCondition(BoundaryConditionType::Neumann, 0.0);
BoundaryCondition rightBoundaryCondition(BoundaryConditionType::Neumann, 0.0);
BoundaryCondition1D boundaryConditions(leftBoundaryCondition, rightBoundaryCondition);
pde::GpuDoublePdeInputData1D data(initialCondition, grid, velocity, diffusion, dt, solverType, SpaceDiscretizerType::Centered, boundaryConditions);
pde::dwave1D solver(data);
solver.Advance(steps);
const auto solution = solver.solution->columns[0]->Get();I wrote a simple python script for plotting the results:
- Numerical instability of the centered difference scheme (regardless of the time solver)
- Numerical diffusion induced by the Upwind scheme and solved by the Lax-Wendroff scheme
- Numerical instability of Explicit Euler scheme
- Numerical instability of Gauss-Legendre (4th order diagonally implicit Runge-Kutta)
- Numerical instability of Implicit/Explicit Euler scheme
cl::dvec xGrid = cl::LinSpace<MemorySpace::Device, MathDomain::Double>(0.0f, 1.0f, 32u);
cl::dvec yGrid = cl::LinSpace<MemorySpace::Device, MathDomain::Double>(0.0f, 1.0f, 32u);
double dt = 1e-5;
double xVelocity = .02;
double yVelocity = .05;
double diffusion = 1.0;
auto _xGrid = xGrid.Get();
auto _yGrid = yGrid.Get();
std::vector<double> _initialCondition(xGrid.size() * yGrid.size());
for (unsigned j = 0; j < _yGrid.size(); ++j)
for (unsigned i = 0; i < _xGrid.size(); ++i)
_initialCondition[i + _xGrid.size() * j] = exp(-_xGrid[i] * _xGrid[i] - _yGrid[i] * _yGrid[i]);
cl::dmat initialCondition(_initialCondition, xGrid.size(), yGrid.size());
pde::GpuDoublePdeInputData2D data(initialCondition, xGrid, yGrid, xVelocity, yVelocity, diffusion, dt, solverType, SpaceDiscretizerType::Centered, boundaryConditions);
pde::dad2D solver(data);
const auto solution = solver.solution->columns[0]->Get(); cl::dvec xGrid = cl::LinSpace<MemorySpace::Device, MathDomain::Double>(0.0f, 1.0f, 32u);
cl::dvec yGrid = cl::LinSpace<MemorySpace::Device, MathDomain::Double>(0.0f, 1.0f, 32u);
double dt = 1e-5;
double xVelocity = .02;
auto _xGrid = xGrid.Get();
auto _yGrid = yGrid.Get();
std::vector<double> _initialCondition(xGrid.size() * yGrid.size());
for (unsigned j = 0; j < _yGrid.size(); ++j)
for (unsigned i = 0; i < _xGrid.size(); ++i)
_initialCondition[i + _xGrid.size() * j] = exp(-_xGrid[i] * _xGrid[i] - _yGrid[i] * _yGrid[i]);
cl::dmat initialCondition(_initialCondition, xGrid.size(), yGrid.size());
pde::GpuDoublePdeInputData2D data(initialCondition, xGrid, yGrid, xVelocity, 0.0, 0.0, dt, solverType, SpaceDiscretizerType::Centered, boundaryConditions);
pde::dwave2D solver(data);
const auto solution = solver.solution->columns[0]->Get();- Lax-Wendroff
- Crank-Nicolson
- Explicit Euler
