Skip to content

GitLab

O
openfpm_numerics
Commit 7f8f9cf2 authored by incardon's avatar incardon
Browse files
Options

Adding CMakeList.txt

parent f0f5fc63
Hide whitespace changes
Inline Side-by-side
Showing with 458 additions and 324 deletions
CMakeLists.txt 0 → 100644
View file @ 7f8f9cf2
cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
project(openfpm_numerics LANGUAGES C CXX CUDA)
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/CMakeFiles/)
set(BOOST_INCLUDE ${Boost_INCLUDE_DIR} CACHE PATH "Include directory for BOOST")
set(PETSC_ROOT CACHE PATH "If compiling with linear algebra indicate the PETSC root directory")
set(PARMETIS_ROOT CACHE PATH "Parmetis root directory")
set(METIS_ROOT CACHE PATH "Metis root directory")
set(LIBHILBERT_ROOT CACHE PATH "LibHilbert root path")
set(HDF5_ROOT CACHE PATH "HDF5 root path")
set(EIGEN3_ROOT CACHE PATH "Eigen3 include path")
set(LIBHILBERT_ROOT CACHE PATH "LibHilbert root path")
set(SE_CLASS1 CACHE BOOL "Activate compilation with SE_CLASS1")
set(SE_CLASS2 CACHE BOOL "Activate compilation with SE_CLASS2")
set(SE_CLASS3 CACHE BOOL "Activate compilation with SE_CLASS3")
set(PROFILE_WITH_SCOREP CACHE BOOL "Enable profiling with scorep")
set(ENV{PETSC_DIR} ${PETSC_ROOT})
set(ENV{HDF5_ROOT} ${HDF5_ROOT})
set(ENV{EIGEN3_ROOT} ${EIGEN3_ROOT})
set(METIS_DIR ${METIS_ROOT})
set(PARMETIS_DIR ${PARMETIS_ROOT})
set(ENV{PATH} "$ENV{PATH}:${HDF5_ROOT}/bin")
set(HDF5_PREFER_PARALLEL TRUE)
find_package(Boost 1.52.0 REQUIRED unit_test_framework iostreams program_options)
find_package(CUDA)
find_package(MPI REQUIRED)
find_package(PETSc)
find_package(HDF5 REQUIRED)
find_package(Eigen3)
find_package(LibHilbert REQUIRED)
find_package(Metis REQUIRED)
find_package(ParMetis REQUIRED)
if(PROFILE_WITH_SCOREP)
set(CMAKE_CXX_COMPILER_LAUNCHER "scorep")
set(CMAKE_CC_COMPILER_LAUNCHER "scorep")
set(CMAKE_CUDA_COMPILER_LAUNCHER "scorep")
endif()
if(CUDA_FOUND)
set(OPENFPM_INIT_FILE "initialize/initialize_wrapper_cuda.cu")
else()
set(OPENFPM_INIT_FILE "initialize/initialize_wrapper_cuda.cpp")
endif()
###### CONFIG.h FILE ######
if(SE_CLASS1)
set(DEFINE_SE_CLASS1 "#define SE_CLASS1")
endif()
if(SE_CLASS2)
set(DEFINE_SE_CLASS2 "#define SE_CLASS2")
endif()
if(SE_CLASS3)
set(DEFINE_SE_CLASS3 "#define SE_CLASS3")
endif()
if(CUDA_FOUND)
set(DEFINE_CUDA_GPU "#define CUDA_GPU")
endif()
if (METIS_FOUND)
set(DEFINE_HAVE_METIS "#define HAVE_METIS 1")
else()
message( FATAL_ERROR "Metis is required in order to install OpenFPM" )
endif()
if (PARMETIS_FOUND)
set(DEFINE_HAVE_PARMETIS "#define HAVE_PARMETIS 1")
else()
message( FATAL_ERROR "ParMetis is required in order to install OpenFPM")
endif()
if(MPI_FOUND)
set(DEFINE_HAVE_MPI "#define HAVE_MPI")
else()
message( FATAL_ERROR "MPI is required in order to install OpenFPM" )
endif()
if (Boost_FOUND)
set(DEFINE_HAVE_BOOST "#define HAVE_BOOST")
set(DEFINE_HAVE_BOOST_IOSTREAMS "#define HAVE_BOOST_IOSTREAMS")
set(DEFINE_HAVE_BOOST_PROGRAM_OPTIONS "#define HAVE_BOOST_PROGRAM_OPTIONS")
set(DEFINE_HAVE_BOOST_UNIT_TEST_FRAMEWORK "#define HAVE_BOOST_UNIT_TEST_FRAMEWORK")
else()
message( FATAL_ERROR "BOOST is required in order to install OpenFPM" )
endif()
if(HDF5_FOUND)
if (HDF5_IS_PARALLEL)
set(DEFINE_HAVE_HDF5 "#define HAVE_HDF5")
else()
message( STATUS "HDF5 found ${HDF5_INCLUDE_DIRS} does not have parallel support, OpenFPM require it" )
endif()
else()
message( FATAL_ERROR "HDF5 with parallel support is required in order to install OpenFPM" )
endif()
if(EIGEN_FOUND)
set(DEFINE_HAVE_EIGEN "#define HAVE_EIGEN")
endif()
if(LIBHILBERT_FOUND)
set(DEFINE_HAVE_LIBHILBERT "#define HAVE_LIBHILBERT 1")
else()
message( FATAL_ERROR "LibHilbert is required in order to install OpenFPM")
endif()
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/src/config/config_cmake.h.in ${CMAKE_CURRENT_SOURCE_DIR}/src/config/config.h)
include_directories(SYSTEM ${MPI_INCLUDE_PATH})
add_subdirectory (src)
configure.ac
View file @ 7f8f9cf2
...@@ -161,7 +161,7 @@ have_quad_head=no ...@@ -161,7 +161,7 @@ have_quad_head=no
AC_CHECK_LIB(quadmath, sinq, [have_quad_lib=yes], []) AC_CHECK_LIB(quadmath, sinq, [have_quad_lib=yes], [])
AC_CHECK_HEADER(quadmath.h,[have_quad_head=yes],[]) AC_CHECK_HEADER(quadmath.h,[have_quad_head=yes],[])
if [x"have_quad_math" == x"yes" $&& x"have_quad_math" == x"yes" ]; then if test x"have_quad_math" == x"yes"; then
AC_DEFINE(HAVE_LIBQUADMATH,[],[Have quad math lib]) AC_DEFINE(HAVE_LIBQUADMATH,[],[Have quad math lib])
LIBQUADMATH=" -lquadmath " LIBQUADMATH=" -lquadmath "
fi fi
......
src/CMakeLists.txt 0 → 100644
View file @ 7f8f9cf2
cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
########################### Executables
add_executable(numerics main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp)
###########################
target_compile_options(numerics PUBLIC $<$<COMPILE_LANGUAGE:CUDA>:-Xcudafe "--display_error_number --diag_suppress=2885 --diag_suppress=2887 --diag_suppress=2888 --diag_suppress=186 --diag_suppress=111" --expt-extended-lambda>)
target_include_directories (numerics PUBLIC ${CUDA_INCLUDE_DIRS})
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../../openfpm_devices/src/)
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../../openfpm_vcluster/src/)
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../../src/)
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../../openfpm_data/src/)
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../../openfpm_io/src/)
target_include_directories (numerics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/config)
target_include_directories (numerics PUBLIC ${PETSC_INCLUDES})
target_include_directories (numerics PUBLIC ${PARMETIS_ROOT}/include)
target_include_directories (numerics PUBLIC ${METIS_ROOT}/include)
target_include_directories (numerics PUBLIC ${HDF5_ROOT}/include)
target_include_directories (numerics PUBLIC ${LIBHILBERT_INCLUDE_DIRS})
target_link_libraries(numerics ${Boost_LIBRARIES})
target_link_libraries(numerics -L${METIS_ROOT}/lib metis)
target_link_libraries(numerics -L${PARMETIS_ROOT}/lib parmetis)
target_link_libraries(numerics -L${HDF5_ROOT}/lib hdf5 hdf5_hl)
target_link_libraries(numerics -L${LIBHILBERT_LIBRARY_DIRS} ${LIBHILBERT_LIBRARIES})
if(PETSC_FOUND)
target_link_libraries(numerics -L${PETSC_LIB_DIR} ${PETSC_LIBRARIES})
endif()
# Request that particles be built with -std=c++11
# As this is a public compile feature anything that links to particles
# will also build with -std=c++11
target_compile_features(numerics PUBLIC cxx_std_11)
target_link_libraries(numerics ${MPI_C_LIBRARIES})
install(FILES Matrix/SparseMatrix.hpp
Matrix/SparseMatrix_Eigen.hpp
Matrix/SparseMatrix_petsc.hpp
DESTINATION openfpm_numerics/include/Matrix)
install(FILES Vector/Vector_eigen.hpp
Vector/Vector_petsc.hpp
Vector/Vector_util.hpp
Vector/Vector.hpp
DESTINATION openfpm_numerics/Vector )
install(FILES Solvers/umfpack_solver.hpp Solvers/petsc_solver.hpp
DESTINATION openfpm_numerics/include/Solvers )
install(FILES util/petsc_util.hpp
util/linalgebra_lib.hpp
util/util_num.hpp
util/grid_dist_testing.hpp
DESTINATION openfpm_numerics/include/util )
install(FILES FiniteDifference/Average.hpp
FiniteDifference/Derivative.hpp
FiniteDifference/FD_util_include.hpp
FiniteDifference/eq.hpp
FiniteDifference/FDScheme.hpp
FiniteDifference/Laplacian.hpp
FiniteDifference/mul.hpp
FiniteDifference/sum.hpp
DESTINATION openfpm_numerics/include/FiniteDifference )
install(FILES FiniteDifference/util/common.hpp
DESTINATION openfpm_numerics/include/FiniteDifference/util )
install(FILES PSE/Kernels.hpp PSE/Kernels_test_util.hpp
DESTINATION openfpm_numerics/include/PSE )
install(FILES Operators/Vector/vector_dist_operators_extensions.hpp
Operators/Vector/vector_dist_operators.hpp
Operators/Vector/vector_dist_operators_apply_kernel.hpp
Operators/Vector/vector_dist_operators_functions.hpp
Operators/Vector/vector_dist_operator_assign.hpp
DESTINATION openfpm_numerics/include/Operators/Vector )
install(FILES Draw/DrawParticles.hpp
Draw/PointIterator.hpp
Draw/PointIteratorSkin.hpp
DESTINATION openfpm_numerics/include/Draw )
install(FILES interpolation/interpolation.hpp
interpolation/mp4_kernel.hpp
interpolation/z_spline.hpp
DESTINATION openfpm_numerics/include/interpolation )
#if(BUILD_TESTING)
# add_executable(particle_test test.cu)
# set_target_properties(particle_test PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
# target_link_libraries(particle_test PRIVATE particles)
# add_test(NAME particles_10k COMMAND particle_test 10000 )
# add_test(NAME particles_256k COMMAND particle_test 256000 )
# if(APPLE)
# We need to add the default path to the driver (libcuda.dylib) as an rpath,
# so that the static cuda runtime can find it at runtime.
# set_property(TARGET particle_test PROPERTY BUILD_RPATH ${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES})
# endif()
#endif()
src/Draw/DrawParticles.hpp
View file @ 7f8f9cf2
...@@ -43,8 +43,9 @@ public: ...@@ -43,8 +43,9 @@ public:
* \return an iterator to the selected particles * \return an iterator to the selected particles
* *
*/ */
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition> template<unsigned int dim, typename T, typename vd_type>
DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd, static PointIteratorSkin<dim,T,typename vd_type::Decomposition_type>
DrawSkin(vd_type & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
Box<dim,T> & domain, Box<dim,T> & domain,
Box<dim,T> & sub_A, Box<dim,T> & sub_A,
...@@ -55,7 +56,7 @@ public: ...@@ -55,7 +56,7 @@ public:
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
bc[i] = NON_PERIODIC; bc[i] = NON_PERIODIC;
return PointIteratorSkin<dim,T,Decomposition>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A, sub_B, bc); return PointIteratorSkin<dim,T,typename vd_type::Decomposition_type>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A, sub_B, bc);
} }
...@@ -79,8 +80,9 @@ public: ...@@ -79,8 +80,9 @@ public:
* \return an iterator to the selected particles * \return an iterator to the selected particles
* *
*/ */
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIteratorSkin<dim,T,Decomposition> template<unsigned int dim, typename T, typename vd_type>
DrawSkin(vector_dist<dim,T,aggr,Decomposition> & vd, static PointIteratorSkin<dim,T,typename vd_type::Decomposition_type>
DrawSkin(vd_type & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
Box<dim,T> & domain, Box<dim,T> & domain,
openfpm::vector<Box<dim,T>> & sub_A, openfpm::vector<Box<dim,T>> & sub_A,
...@@ -91,7 +93,7 @@ public: ...@@ -91,7 +93,7 @@ public:
for (size_t i = 0 ; i < dim ; i++) for (size_t i = 0 ; i < dim ; i++)
bc[i] = NON_PERIODIC; bc[i] = NON_PERIODIC;
PointIteratorSkin<dim,T,Decomposition> it(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A.get(0), sub_B, bc); PointIteratorSkin<dim,T,typename vd_type::Decomposition_type> it(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub_A.get(0), sub_B, bc);
for (size_t i = 1 ; i < sub_A.size() ; i++) for (size_t i = 1 ; i < sub_A.size() ; i++)
it.addBoxA(Box<dim,T>(sub_A.get(i))); it.addBoxA(Box<dim,T>(sub_A.get(i)));
...@@ -117,13 +119,13 @@ public: ...@@ -117,13 +119,13 @@ public:
* \return an iterator to the selected particles * \return an iterator to the selected particles
* *
*/ */
template<unsigned int dim, typename T, typename aggr, typename Decomposition> static PointIterator<dim,T,Decomposition> template<unsigned int dim, typename T, typename vd_type> static PointIterator<dim,T,typename vd_type::Decomposition_type>
DrawBox(vector_dist<dim,T,aggr,Decomposition> & vd, DrawBox(vd_type & vd,
size_t (& sz)[dim], size_t (& sz)[dim],
Box<dim,T> & domain, Box<dim,T> & domain,
Box<dim,T> & sub) Box<dim,T> & sub)
{ {
return PointIterator<dim,T,Decomposition>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub); return PointIterator<dim,T,typename vd_type::Decomposition_type>(vd.getDecomposition(),sz,vd.getDecomposition().getDomain(),sub);
} }
}; };
......
src/Draw/DrawParticles_unit_tests.hpp
View file @ 7f8f9cf2
...@@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE(point_iterator) ...@@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE(point_iterator)
//! [DrawBox_example] //! [DrawBox_example]
Vcluster & v_cl = create_vcluster(); Vcluster<> & v_cl = create_vcluster();
v_cl.sum(cnt); v_cl.sum(cnt);
v_cl.execute(); v_cl.execute();
...@@ -115,7 +115,7 @@ BOOST_AUTO_TEST_CASE(point_iterator_skin) ...@@ -115,7 +115,7 @@ BOOST_AUTO_TEST_CASE(point_iterator_skin)
++p; ++p;
} }
Vcluster & v_cl = create_vcluster(); Vcluster<> & v_cl = create_vcluster();
v_cl.sum(cnt); v_cl.sum(cnt);
v_cl.execute(); v_cl.execute();
......
src/Equations/stoke_flow_eq_3d.hpp
View file @ 7f8f9cf2
...@@ -8,6 +8,10 @@ ...@@ -8,6 +8,10 @@
#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_STOKE_FLOW_EQ_3D_HPP_ #ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_STOKE_FLOW_EQ_3D_HPP_
#define OPENFPM_NUMERICS_SRC_EQUATIONS_STOKE_FLOW_EQ_3D_HPP_ #define OPENFPM_NUMERICS_SRC_EQUATIONS_STOKE_FLOW_EQ_3D_HPP_
constexpr unsigned int x = 0;
constexpr unsigned int y = 1;
constexpr unsigned int z = 2;
// Model the equations // Model the equations
constexpr unsigned int v[] = {0,1,2}; constexpr unsigned int v[] = {0,1,2};
......
src/FiniteDifference/FDScheme.hpp
View file @ 7f8f9cf2
...@@ -310,7 +310,7 @@ private: ...@@ -310,7 +310,7 @@ private:
// Indicate all the non zero colums // Indicate all the non zero colums
// This check can be done only on single processor // This check can be done only on single processor
Vcluster & v_cl = create_vcluster(); Vcluster<> & v_cl = create_vcluster();
if (v_cl.getProcessingUnits() == 1) if (v_cl.getProcessingUnits() == 1)
{ {
openfpm::vector<unsigned> nz_cols; openfpm::vector<unsigned> nz_cols;
...@@ -604,7 +604,7 @@ private: ...@@ -604,7 +604,7 @@ private:
*/ */
void construct_gmap() void construct_gmap()
{ {
Vcluster & v_cl = create_vcluster(); Vcluster<> & v_cl = create_vcluster();
// Calculate the size of the local domain // Calculate the size of the local domain
size_t sz = g_map.getLocalDomainSize(); size_t sz = g_map.getLocalDomainSize();
......
src/FiniteDifference/FDScheme_unit_tests.hpp → src/FiniteDifference/FDScheme_unit_tests.cpp
View file @ 7f8f9cf2
...@@ -8,6 +8,9 @@ ...@@ -8,6 +8,9 @@
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_UNIT_TESTS_HPP_ #ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_UNIT_TESTS_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_UNIT_TESTS_HPP_ #define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_UNIT_TESTS_HPP_
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#include "FiniteDifference/Derivative.hpp" #include "FiniteDifference/Derivative.hpp"
#include "FiniteDifference/Laplacian.hpp" #include "FiniteDifference/Laplacian.hpp"
#include "Decomposition/CartDecomposition.hpp" #include "Decomposition/CartDecomposition.hpp"
......
src/FiniteDifference/eq.hpp
View file @ 7f8f9cf2
...@@ -15,6 +15,7 @@ ...@@ -15,6 +15,7 @@
//#define NON_PERIODIC false //#define NON_PERIODIC false
#include "util/util_num.hpp" #include "util/util_num.hpp"
#include "Matrix/SparseMatrix.hpp"
/*! \brief Equation /*! \brief Equation
* *
......
src/FiniteDifference/eq_unit_test.hpp deleted 100644 → 0
View file @ f0f5fc63
/*
* eq_unit_test.hpp
*
* Created on: Oct 13, 2015
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_EQ_UNIT_TEST_HPP_
#include "Laplacian.hpp"
#include "FiniteDifference/eq.hpp"
#include "FiniteDifference/sum.hpp"
#include "FiniteDifference/mul.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "data_type/aggregate.hpp"
#include "FiniteDifference/FDScheme.hpp"
BOOST_AUTO_TEST_SUITE( eq_test_suite )
//! [Definition of the system]
struct lid_nn
{
// dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system v_x, v_y, P so a total of 3
static const unsigned int nvar = 3;
// boundary conditions PERIODIC OR NON_PERIODIC
static const bool boundary[];
// type of space float, double, ...
typedef float stype;
// type of base grid, it is the distributed grid that will store the result
// Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack it is the only possible choice
typedef SparseMatrix<double,int> SparseMatrix_type;
// type of Vector for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack it is the only possible choice
typedef Vector<double> Vector_type;
// Define that the underline grid where we discretize the system of equation is staggered
static const int grid_type = STAGGERED_GRID;
};
const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
//! [Definition of the system]
//! [Definition of the equation of the system in the bulk and at the boundary]
// Constant Field
struct eta
{
typedef void const_field;
static float val() {return 1.0;}
};
// Convenient constants
constexpr unsigned int v[] = {0,1};
constexpr unsigned int P = 2;
constexpr unsigned int ic = 2;
// Create field that we have v_x, v_y, P
typedef Field<v[x],lid_nn> v_x;
typedef Field<v[y],lid_nn> v_y;
typedef Field<P,lid_nn> Prs;
// Eq1 V_x
typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx;
typedef D<x,Prs,lid_nn> p_x;
typedef minus<p_x,lid_nn> m_p_x;
typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq;
// Eq2 V_y
typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
typedef D<y,Prs,lid_nn> p_y;
typedef minus<p_y,lid_nn> m_p_y;
typedef sum<eta_lap_vy,m_p_y,lid_nn> vy_eq;
// Eq3 Incompressibility
typedef D<x,v_x,lid_nn,FORWARD> dx_vx;
typedef D<y,v_y,lid_nn,FORWARD> dy_vy;
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
// Equation for boundary conditions
/* Consider the staggered cell
*
\verbatim
+--$--+
| |
# * #
| |
0--$--+
# = velocity(y)
$ = velocity(x)
* = pressure
\endverbatim
*
*
* If we want to impose v_y = 0 on 0 we have to interpolate between # of this cell
* and # of the previous cell on y, (Average) or Avg operator
*
*/
// Directional Avg
typedef Avg<x,v_y,lid_nn> avg_vy;
typedef Avg<y,v_x,lid_nn> avg_vx;
typedef Avg<x,v_y,lid_nn,FORWARD> avg_vy_f;
typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
//! [Definition of the equation of the system in the bulk and at the boundary]
template<typename solver_type,typename lid_nn> void lid_driven_cavity_2d()
{
Vcluster & v_cl = create_vcluster();
if (v_cl.getProcessingUnits() > 3)
return;
//! [lid-driven cavity 2D]
// velocity in the grid is the property 0, pressure is the property 1
constexpr int velocity = 0;
constexpr int pressure = 1;
// Domain, a rectangle
Box<2,float> domain({0.0,0.0},{3.0,1.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=256 and y=64
long int sz[] = {256,64};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
// We need one more point on the left and down part of the domain
// This is given by the boundary conditions that we impose, the
// reason is mathematical in order to have a well defined system
// and cannot be discussed here
Padding<2> pd({1,1},{0,0});
// Distributed grid that store the solution
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
// It is the maximum extension of the stencil
Ghost<2,long int> stencil_max(1);
// Finite difference scheme
FDScheme<lid_nn> fd(pd, stencil_max, domain,g_dist);
// Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
// exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
// mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
// solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
// we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
// Here we impose the Eq1 and Eq2
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
// v_x and v_y
// Imposing B1
fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
// Imposing B2
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
// Imposing B3
fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
// Imposing B4
fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
// When we pad the grid, there are points of the grid that are not
// touched by the previous condition. Mathematically this lead