diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 61c88a9c010e6d500e89acb84c9328339d69c67c..10aca02ef7d30d726a58df0e903ca7aebcd821bf 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -200,4 +200,159 @@ else()
cube_1e3_neumann_pcs_0_ts_1_t_1.000000_1.vtu pressure pressure
cube_1e3_neumann_pcs_0_ts_1_t_1.000000_2.vtu pressure pressure
)
+
+ # Single core
+ # CUBE 1x1x1 GROUNDWATER FLOW TESTS
+ foreach(mesh_size 1e0 1e1 1e2 1e3)
+ AddTest(
+ NAME GroundWaterFlowProcess_cube_1x1x1_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE_ARGS cube_${mesh_size}.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA cube_${mesh_size}_pcs_0_ts_1_t_1_0.vtu Linear_1_to_minus1 pressure
+ )
+
+ AddTest(
+ NAME GroundWaterFlowProcess_cube_1x1x1_Neumann_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-1 RELTOL 1e-1
+ DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1_t_1_0.vtu D1_left_front_N1_right pressure
+ )
+ endforeach()
+
+
+ foreach(mesh_size 1e4 2e4 3e4 4e4 5e4 1e5 1e6)
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_cube_1x1x1_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE_ARGS cube_${mesh_size}.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-7 RELTOL 1e-7
+ DIFF_DATA cube_${mesh_size}_pcs_0_ts_1_t_1_0.vtu Linear_1_to_minus1 pressure
+ )
+
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_cube_1x1x1_Neumann_${mesh_size}
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE_ARGS cube_${mesh_size}_neumann.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-2 RELTOL 1e-2
+ DIFF_DATA cube_${mesh_size}_neumann_pcs_0_ts_1_t_1_0.vtu D1_left_front_N1_right pressure
+ )
+ endforeach()
+
+ # SQUARE 1x1 GROUNDWATER FLOW TESTS
+ foreach(mesh_size 1e0 1e1 1e2 1e3 1e4)
+ AddTest(
+ NAME GroundWaterFlowProcess_square_1x1_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE_ARGS square_${mesh_size}.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-13 RELTOL 1e-13
+ DIFF_DATA square_${mesh_size}_pcs_0_ts_1_t_1_0.vtu Linear_1_to_minus1 pressure
+ )
+
+ AddTest(
+ NAME GroundWaterFlowProcess_square_1x1_Neumann_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-1 RELTOL 1e-1
+ DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1_t_1_0.vtu D1_left_bottom_N1_right pressure
+ )
+ endforeach()
+
+ foreach(mesh_size 1e5 1e6)
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_square_1x1_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE_ARGS square_${mesh_size}.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-7 RELTOL 1e-7
+ DIFF_DATA square_${mesh_size}_pcs_0_ts_1_t_1_0.vtu Linear_1_to_minus1 pressure
+ )
+
+ AddTest(
+ NAME LARGE_GroundWaterFlowProcess_square_1x1_Neumann_${mesh_size}
+ PATH Elliptic/square_1x1_GroundWaterFlow
+ EXECUTABLE_ARGS square_${mesh_size}_neumann.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-02 RELTOL 1e-02
+ DIFF_DATA square_${mesh_size}_neumann_pcs_0_ts_1_t_1_0.vtu D1_left_bottom_N1_right pressure
+ )
+ endforeach()
+
+ # LINE 1 GROUNDWATER FLOW TESTS
+ foreach(mesh_size 1e1)
+ AddTest(
+ NAME GroundWaterFlowProcess_line_1_${mesh_size}
+ PATH Elliptic/line_1_GroundWaterFlow
+ EXECUTABLE_ARGS line_${mesh_size}.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA line_${mesh_size}_pcs_0_ts_1_t_1_0.vtu Linear_1_to_minus1 pressure
+ )
+
+ AddTest(
+ NAME GroundWaterFlowProcess_line_1_Neumann_${mesh_size}
+ PATH Elliptic/line_1_GroundWaterFlow
+ EXECUTABLE_ARGS line_${mesh_size}_neumann.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-14 RELTOL 1e-14
+ DIFF_DATA line_${mesh_size}_neumann_pcs_0_ts_1_t_1_0.vtu D1_left_N1_right pressure
+ )
+ endforeach()
+
+ AddTest(
+ NAME TES_zeolite_discharge_small
+ PATH Parabolic/TES/1D
+ EXECUTABLE_ARGS tes-1D-zeolite-discharge-small.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-16 RELTOL 1e-16
+ DIFF_DATA
+ tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000_0.vtu pressure pressure
+ tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000_0.vtu temperature temperature
+ tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000_0.vtu v_mass_frac v_mass_frac
+# tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000_0.vtu solid_density solid_density
+ )
+
+ AddTest(
+ NAME LARGE_TES_zeolite_discharge
+ PATH Parabolic/TES/1D
+ EXECUTABLE_ARGS tes-1D-zeolite-discharge-large.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-16 RELTOL 1e-16
+ DIFF_DATA
+ tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000_0.vtu pressure pressure
+ tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000_0.vtu temperature temperature
+ tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000_0.vtu v_mass_frac v_mass_frac
+# tes_zeolite_discharge_large_pcs_0_ts_28_t_1_0.vtu solid_density solid_density
+ )
endif()
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
index 59d0ba603712f55c16e823417c2e1498cef6e3c9..961678f85bf8d423d9cf039d5cb2ec1463f0a5be 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -107,6 +107,12 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
}
else if(reason == KSP_DIVERGED_ITS)
{
+ const char *ksp_type;
+ const char *pc_type;
+ KSPGetType(_solver, &ksp_type);
+ PCGetType(_pc, &pc_type);
+ PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver %s with %s preconditioner",
+ ksp_type, pc_type);
PetscPrintf(PETSC_COMM_WORLD, "\nWarning: maximum number of iterations reached.\n");
}
else
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index 2695a1b755a75aca736e6fae3449d82e5deaa351..6d380a12cee3ab88047f48690b1c75225d9d7e99 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -91,6 +91,10 @@ void PETScMatrix::setRowsColumnsZero(std::vector<PetscInt> const& row_pos)
// and thus improves performance for very large process counts.
// See PETSc doc about MAT_NO_OFF_PROC_ZERO_ROWS.
MatSetOption(_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
+
+ // Keep the non-zero pattern for the assignment operator.
+ MatSetOption(_A, MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
+
if(nrows>0)
MatZeroRows(_A, nrows, &row_pos[0], one, PETSC_NULL, PETSC_NULL);
else
diff --git a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
index c48980a5e1db4ae6eaa7064b69a3d7a2835d9f16..f578847886d9860d67a668de05ca341f486c6d3b 100644
--- a/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
+++ b/MeshLib/IO/MPI_IO/NodePartitionedMeshReader.cpp
@@ -418,7 +418,6 @@ NodePartitionedMeshReader::newMesh(
mesh_name + std::to_string(_mpi_comm_size),
mesh_nodes, glb_node_ids, mesh_elems,
MeshLib::Properties(),
- _mesh_info.regular_elements,
_mesh_info.global_base_nodes,
_mesh_info.global_nodes,
_mesh_info.base_nodes,
diff --git a/MeshLib/IO/readMeshFromFile.cpp b/MeshLib/IO/readMeshFromFile.cpp
index f76cb53320e6fbc26372598dfd084ffc94490f57..9f0fcf9f7102e4a9168eea2909ae649fded3808c 100644
--- a/MeshLib/IO/readMeshFromFile.cpp
+++ b/MeshLib/IO/readMeshFromFile.cpp
@@ -45,10 +45,30 @@ namespace IO
MeshLib::Mesh* readMeshFromFile(const std::string &file_name)
{
#ifdef USE_PETSC
- MeshLib::IO::NodePartitionedMeshReader read_pmesh(PETSC_COMM_WORLD);
- const std::string file_name_base = BaseLib::dropFileExtension(file_name);
- return read_pmesh.read(file_name_base);
+ int world_size;
+ MPI_Comm_size(MPI_COMM_WORLD, &world_size);
+ if (world_size > 1)
+ {
+ MeshLib::IO::NodePartitionedMeshReader read_pmesh(PETSC_COMM_WORLD);
+ const std::string file_name_base = BaseLib::dropFileExtension(file_name);
+ return read_pmesh.read(file_name_base);
+ }
+ else if (world_size == 1)
+ {
+ MeshLib::Mesh* mesh = readMeshFromFileSerial(file_name);
+ MeshLib::NodePartitionedMesh* part_mesh
+ = new MeshLib::NodePartitionedMesh(*mesh);
+ delete mesh;
+ return part_mesh;
+ }
+ return nullptr;
#else
+ return readMeshFromFileSerial(file_name);
+#endif
+}
+
+MeshLib::Mesh* readMeshFromFileSerial(const std::string &file_name)
+{
if (BaseLib::hasFileExtension("msh", file_name))
{
MeshLib::IO::Legacy::MeshIO meshIO;
@@ -60,7 +80,8 @@ MeshLib::Mesh* readMeshFromFile(const std::string &file_name)
ERR("readMeshFromFile(): Unknown mesh file format in file %s.", file_name.c_str());
return nullptr;
-#endif
}
+
+
} // end namespace IO
} // end namespace MeshLib
diff --git a/MeshLib/IO/readMeshFromFile.h b/MeshLib/IO/readMeshFromFile.h
index 7fa9e8fa464ded7ae3a1d0cde86278332db7f497..15801bf3747ce02609458d48dd350eb8f5b7280f 100644
--- a/MeshLib/IO/readMeshFromFile.h
+++ b/MeshLib/IO/readMeshFromFile.h
@@ -27,6 +27,7 @@ class Mesh;
namespace IO
{
+MeshLib::Mesh* readMeshFromFileSerial(const std::string &file_name);
MeshLib::Mesh* readMeshFromFile(const std::string &file_name);
}
}
diff --git a/MeshLib/Node.h b/MeshLib/Node.h
index 410f25614f997524b5e563a0add6a9cd812c7385..8b3b573d8ef05634281295e187d974090205a351 100644
--- a/MeshLib/Node.h
+++ b/MeshLib/Node.h
@@ -38,6 +38,7 @@ class Node final : public MathLib::Point3dWithID
{
/* friend classes: */
friend class Mesh;
+ friend class NodePartitionedMesh;
friend class MeshRevision;
friend class MeshLayerMapper;
friend class ApplicationUtils::NodeWiseMeshPartitioner;
diff --git a/MeshLib/NodePartitionedMesh.h b/MeshLib/NodePartitionedMesh.h
index eab5679910e1931811e190ada2590e77302bbeb4..5cf60f536fb1892c34ef12b963109d2e9c25f343 100644
--- a/MeshLib/NodePartitionedMesh.h
+++ b/MeshLib/NodePartitionedMesh.h
@@ -29,6 +29,36 @@ namespace MeshLib
class NodePartitionedMesh : public Mesh
{
public:
+ // Copy a global mesh for the case of the thread number is one,
+ // i.e the gobal mesh is not partitioned.
+ // \param mesh The gobal mesh
+ explicit NodePartitionedMesh(const Mesh& mesh)
+ : Mesh(mesh), _global_node_ids(mesh.getNumberOfNodes()),
+ _n_global_base_nodes(mesh.getNumberOfBaseNodes()),
+ _n_global_nodes(mesh.getNumberOfNodes()),
+ _n_active_base_nodes(mesh.getNumberOfBaseNodes()),
+ _n_active_nodes(mesh.getNumberOfNodes())
+ {
+ const auto& mesh_nodes = mesh.getNodes();
+ for (std::size_t i = 0; i < _nodes.size(); i++)
+ {
+ _global_node_ids[i] = _nodes[i]->getID();
+
+ // TODO To add copying of the connected nodes (and elements)
+ // in the copy constructor of class Node in order to
+ // drop the following lines.
+ auto node = _nodes[i];
+ // Copy constructor of Mesh does not copy the connected
+ // nodes to node.
+ if (node->_connected_nodes.size() == 0)
+ {
+ std::copy(mesh_nodes[i]->_connected_nodes.begin(),
+ mesh_nodes[i]->_connected_nodes.end(),
+ std::back_inserter(node->_connected_nodes));
+ }
+ }
+ }
+
/*!
\brief Constructor
\param name Name assigned to the mesh.
@@ -39,8 +69,7 @@ class NodePartitionedMesh : public Mesh
\param glb_node_ids Global IDs of nodes of a partition.
\param elements Vector for elements. Ghost elements are stored
after regular (non-ghost) elements.
- \param n_nghost_elem Number of non-ghost elements, or the start ID of
- the entry of ghost element in the element vector.
+ \param properties Mesh property.
\param n_global_base_nodes Number of the base nodes of the global mesh.
\param n_global_nodes Number of all nodes of the global mesh.
\param n_base_nodes Number of the base nodes.
@@ -52,14 +81,13 @@ class NodePartitionedMesh : public Mesh
const std::vector<std::size_t> &glb_node_ids,
const std::vector<Element*> &elements,
Properties properties,
- const std::size_t n_nghost_elem,
const std::size_t n_global_base_nodes,
const std::size_t n_global_nodes,
const std::size_t n_base_nodes,
const std::size_t n_active_base_nodes,
const std::size_t n_active_nodes)
: Mesh(name, nodes, elements, properties, n_base_nodes),
- _global_node_ids(glb_node_ids), _n_nghost_elem(n_nghost_elem),
+ _global_node_ids(glb_node_ids),
_n_global_base_nodes(n_global_base_nodes),
_n_global_nodes(n_global_nodes),
_n_active_base_nodes(n_active_base_nodes),
@@ -114,12 +142,6 @@ class NodePartitionedMesh : public Mesh
return _n_base_nodes + _n_active_nodes - _n_active_base_nodes;
}
- /// Get the number of non-ghost elements, or the start entry ID of ghost elements in element vector.
- std::size_t getNumberOfNonGhostElements() const
- {
- return _n_nghost_elem;
- }
-
// TODO I guess that is a simplified version of computeSparsityPattern()
/// Get the maximum number of connected nodes to node.
std::size_t getMaximumNConnectedNodesToNode() const
@@ -139,9 +161,6 @@ class NodePartitionedMesh : public Mesh
/// Global IDs of nodes of a partition
std::vector<std::size_t> _global_node_ids;
- /// Number of non-ghost elements, or the ID of the start entry of ghost elements in _elements vector.
- std::size_t _n_nghost_elem;
-
/// Number of the nodes of the global mesh linear interpolations.
std::size_t _n_global_base_nodes;
diff --git a/NumLib/DOF/ComputeSparsityPattern.cpp b/NumLib/DOF/ComputeSparsityPattern.cpp
index 46d2ad519d03acca1792ad15e3b469da6417bf16..cf80de61fe44ea3ce04494bd2ea693fbcb05c2bd 100644
--- a/NumLib/DOF/ComputeSparsityPattern.cpp
+++ b/NumLib/DOF/ComputeSparsityPattern.cpp
@@ -16,14 +16,15 @@
#include "MeshLib/NodePartitionedMesh.h"
GlobalSparsityPattern computeSparsityPatternPETSc(
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh)
{
assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(&mesh));
auto const& npmesh =
*static_cast<MeshLib::NodePartitionedMesh const*>(&mesh);
- auto const max_nonzeroes = npmesh.getMaximumNConnectedNodesToNode();
+ auto const max_nonzeroes = dof_table.getNumberOfComponents()
+ * npmesh.getMaximumNConnectedNodesToNode();
// The sparsity pattern is misused here in the sense that it will only
// contain a single value.
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index af0dc0486c8748c0646ecfedc8e0bf9732091326..b52772d98fef9f2fa115121b1e75e06b34cc4fde 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -94,23 +94,28 @@ void checkGlobalMatrixInterfaceMPI(T_MATRIX &m, T_VECTOR &v)
MathLib::finalizeMatrixAssembly(m);
+ // Test basic assignment operator with an empty T_MATRIX._A
+ T_MATRIX m_c = m;
+ // Test basic assignment operator with an initalized T_MATRIX._A
+ m_c = m;
+
// Multiply by a vector
// v = 1.;
set(v, 1.);
const bool deep_copy = false;
T_VECTOR y(v, deep_copy);
- matMult(m, v, y);
+ matMult(m_c, v, y);
ASSERT_EQ(sqrt(3*(3*3 + 7*7)), norm2(y));
// set a value
- m.set(2 * mrank, 2 * mrank, 5.0);
+ m_c.set(2 * mrank, 2 * mrank, 5.0);
MathLib::finalizeMatrixAssembly(m);
// add a value
- m.add(2 * mrank+1, 2 * mrank+1, 5.0);
- MathLib::finalizeMatrixAssembly(m);
+ m_c.add(2 * mrank+1, 2 * mrank+1, 5.0);
+ MathLib::finalizeMatrixAssembly(m_c);
- matMult(m, v, y);
+ matMult(m_c, v, y);
ASSERT_EQ(sqrt((3*7*7 + 3*12*12)), norm2(y));
}