diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp
index 1f804d46c6ece3796c62a8363ad2af6679ed481b..927fe14b792bd66b017202335523c145747e47f2 100644
--- a/MathLib/LinAlg/LinAlg.cpp
+++ b/MathLib/LinAlg/LinAlg.cpp
@@ -22,6 +22,11 @@ namespace MathLib { namespace LinAlg
// Vector
+void setLocalAccessibleVector(PETScVector const& x)
+{
+ x.setLocalAccessibleVector();
+}
+
void set(PETScVector& x, double const a)
{
VecSet(x.getRawVector(), a);
@@ -177,6 +182,10 @@ namespace MathLib { namespace LinAlg
// Vector
+void setLocalAccessibleVector(EigenVector const& /*x*/)
+{
+}
+
void set(EigenVector& x, double const a)
{
x.getRawVector().setConstant(a);
diff --git a/MathLib/LinAlg/LinAlg.h b/MathLib/LinAlg/LinAlg.h
index 06388524bce6be8f241a152211a724c502cb9678..8577adf1d4e6a879830cf9998b29fd5b1de4d570 100644
--- a/MathLib/LinAlg/LinAlg.h
+++ b/MathLib/LinAlg/LinAlg.h
@@ -148,6 +148,10 @@ namespace LinAlg
// Vector
+/// Set local accessible vector in order to get entries.
+/// Call this before call operator[] or get(...) of x.
+void setLocalAccessibleVector(PETScVector const& x);
+
void set(PETScVector& x, double const a);
void copy(PETScVector const& x, PETScVector& y);
@@ -206,6 +210,16 @@ namespace LinAlg
// Vector
+/**
+ Set local accessible vector in order to get entries. Call this
+ before call operator[] or get(...) of x.
+ The function only has computation if DDC is enabled,
+ e.g. parallel computing. Up to now, Eigen vector is not used for
+ global vectors in parallel computing. Therefore this function is
+ empty in the current status.
+*/
+void setLocalAccessibleVector(EigenVector const& x);
+
void set(EigenVector& x, double const a);
void copy(EigenVector const& x, EigenVector& y);
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
index e4fa6ee8052e066854a8e7fe16ed6c1524906a4c..a563e4e47c9f9cdb117f7aed1ba991216c3f9853 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -83,7 +83,7 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x)
KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
#endif
- KSPSolve(_solver, b.getData(), x.getData());
+ KSPSolve(_solver, b.getRawVector(), x.getRawVector());
KSPConvergedReason reason;
KSPGetConvergedReason(_solver, &reason);
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index e5ea7875ce0168589b5a4142b45b61315bf3920e..12e1a2b10e109965edfed5965cb456dddb26dcc3 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -28,10 +28,13 @@ namespace MathLib
{
PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
{
- if( is_global_size ) {
+ if( is_global_size )
+ {
VecCreate(PETSC_COMM_WORLD, &_v);
VecSetSizes(_v, PETSC_DECIDE, vec_size);
- } else {
+ }
+ else
+ {
// Fix size partitioning
// the size can be associated to specific memory allocation of a matrix
VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, &_v);
@@ -43,8 +46,8 @@ PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size)
PETScVector::PETScVector(const PetscInt vec_size,
const std::vector<PetscInt>& ghost_ids,
const bool is_global_size)
- : _size_ghosts{static_cast<PetscInt>(ghost_ids.size())}
- , _has_ghost_id{true}
+ : _size_ghosts {static_cast<PetscInt>(ghost_ids.size())},
+ _has_ghost_id {true}
{
PetscInt nghosts = static_cast<PetscInt>( ghost_ids.size() );
if ( is_global_size )
@@ -61,6 +64,10 @@ PETScVector::PETScVector(const PetscInt vec_size,
}
config();
+
+ for (PetscInt i=0; i<nghosts; i++)
+ _global_ids2local_ids_ghost.emplace(ghost_ids[i],
+ _size_loc + i);
}
PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
@@ -75,14 +82,15 @@ PETScVector::PETScVector(const PETScVector &existing_vec, const bool deep_copy)
}
PETScVector::PETScVector(PETScVector &&other)
- : _v{std::move(other._v)}
- , _v_loc{std::move(other._v_loc)}
- , _start_rank{other._start_rank}
- , _end_rank{other._end_rank}
- , _size{other._size}
- , _size_loc{other._size_loc}
- , _size_ghosts{other._size_ghosts}
- , _has_ghost_id{other._has_ghost_id}
+ : _v {std::move(other._v)},
+ _v_loc {std::move(other._v_loc)},
+ _start_rank {other._start_rank},
+ _end_rank {other._end_rank},
+ _size {other._size},
+ _size_loc {other._size_loc},
+ _size_ghosts {other._size_ghosts},
+ _has_ghost_id {other._has_ghost_id},
+ _global_ids2local_ids_ghost {other._global_ids2local_ids_ghost}
{}
void PETScVector::config()
@@ -104,7 +112,7 @@ void PETScVector::finalizeAssembly()
}
void PETScVector::gatherLocalVectors( PetscScalar local_array[],
- PetscScalar global_array[])
+ PetscScalar global_array[]) const
{
// Collect vectors from processors.
int size_rank;
@@ -130,7 +138,7 @@ void PETScVector::gatherLocalVectors( PetscScalar local_array[],
}
-void PETScVector::getGlobalVector(PetscScalar u[])
+void PETScVector::getGlobalVector(std::vector<PetscScalar>& u) const
{
#ifdef TEST_MEM_PETSC
@@ -138,10 +146,12 @@ void PETScVector::getGlobalVector(PetscScalar u[])
PetscMemoryGetCurrentUsage(&mem1);
#endif
+ assert(u.size() == _size);
+
PetscScalar *xp = nullptr;
VecGetArray(_v, &xp);
- gatherLocalVectors(xp, u);
+ gatherLocalVectors(xp, u.data());
//This following line may be needed late on
// for a communication load balance:
@@ -156,15 +166,77 @@ void PETScVector::getGlobalVector(PetscScalar u[])
#endif
}
-void PETScVector::copyValues(std::vector<double>& u) const
+void PETScVector::setLocalAccessibleVector() const
+{
+ if (_entry_array.empty())
+ {
+ const PetscInt array_size
+ = _global_ids2local_ids_ghost.size() > 0 ?
+ _size_loc + _size_ghosts: _size;
+ _entry_array.resize(array_size);
+ }
+
+ if ( !_global_ids2local_ids_ghost.empty() )
+ {
+ PetscScalar* loc_x = getLocalVector();
+ std::copy_n(loc_x, _size_loc + _size_ghosts,
+ _entry_array.begin());
+ restoreArray(loc_x);
+ }
+ else
+ getGlobalVector(_entry_array);
+}
+
+void PETScVector::copyValues(std::vector<PetscScalar>& u) const
{
assert(u.size() == (std::size_t) (getLocalSize() + getGhostSize()));
- double* loc_x = getLocalVector();
+ PetscScalar* loc_x = getLocalVector();
std::copy_n(loc_x, getLocalSize() + getGhostSize(), u.begin());
restoreArray(loc_x);
}
+PetscScalar PETScVector::get(const PetscInt idx) const
+{
+ if (_global_ids2local_ids_ghost.size() > 0)
+ {
+ return _entry_array[getLocalIndex(idx)];
+ }
+
+ // Ghost entries, and its original index is 0.
+ const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx);
+ return _entry_array[id_p];
+}
+
+
+std::vector<PetscScalar> PETScVector::get(std::vector<IndexType> const& indices) const
+{
+ std::vector<PetscScalar> local_x(indices.size());
+ // If VecGetValues can get values from different processors,
+ // use VecGetValues(_v, indices.size(), indices.data(),
+ // local_x.data());
+
+ if (_global_ids2local_ids_ghost.size() > 0)
+ {
+ for (std::size_t i=0; i<indices.size(); i++)
+ {
+ local_x[i] = _entry_array[getLocalIndex(indices[i])];
+ }
+ }
+ else
+ {
+ for (std::size_t i=0; i<indices.size(); i++)
+ {
+ // Ghost entries, and its original index is 0.
+ const IndexType id_p = (indices[i] == -_size)
+ ? 0 : std::abs(indices[i]);
+ local_x[i] = _entry_array[id_p];
+ }
+ }
+
+ return local_x;
+}
+
PetscScalar* PETScVector::getLocalVector() const
{
PetscScalar *loc_array;
@@ -176,10 +248,23 @@ PetscScalar* PETScVector::getLocalVector() const
VecGetArray(_v_loc, &loc_array);
}
else
- VecGetArray(_v, &loc_array);
+ VecGetArray(_v, &loc_array);
return loc_array;
}
+PetscInt PETScVector::getLocalIndex(const PetscInt global_index) const
+{
+ if (global_index >= 0) // non-ghost entry.
+ return global_index - _start_rank;
+
+ // A special case for a ghost location with global index equal to
+ // the size of the local vector:
+ PetscInt real_global_index =
+ (-global_index == _size) ? 0 : -global_index;
+
+ return _global_ids2local_ids_ghost.at(real_global_index);
+}
+
void PETScVector::restoreArray(PetscScalar* array) const
{
if (_has_ghost_id)
@@ -221,6 +306,7 @@ void PETScVector::shallowCopy(const PETScVector &v)
_size_loc = v._size_loc;
_size_ghosts = v._size_ghosts;
_has_ghost_id = v._has_ghost_id;
+ _global_ids2local_ids_ghost = v._global_ids2local_ids_ghost;
VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
}
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index cba4692313e73cdc099c0d9cfa1e2048e182c27a..bb5193109d521e3a668c5c0339ca03f7e39ab9c0 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -16,8 +16,10 @@
#pragma once
+#include <map>
#include <string>
#include <vector>
+
#include <petscvec.h>
namespace MathLib
@@ -154,59 +156,41 @@ class PETScVector
VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES);
}
- //! Get several entries
- std::vector<double> get(std::vector<IndexType> const& indices) const
- {
- std::vector<double> local_x(indices.size());
+ // TODO preliminary
+ void setZero() { VecSet(_v, 0.0); }
- VecGetValues(_v, indices.size(), indices.data(), local_x.data());
+ /// Disallow moving.
+ /// \todo This operator should be implemented properly when doing a
+ /// general cleanup of all matrix and vector classes.
+ PETScVector& operator = (PETScVector &&) = delete;
- return local_x;
- }
+ /// Set local accessible vector in order to get entries.
+ /// Call this before call operator[] or get(...).
+ void setLocalAccessibleVector() const;
- // TODO preliminary
- double operator[] (PetscInt idx) const
+ /// Get several entries. setLocalAccessibleVector() must be
+ /// called beforehand.
+ std::vector<PetscScalar> get(std::vector<IndexType> const& indices) const;
+
+ /// Get the value of an entry by [] operator.
+ /// setLocalAccessibleVector() must be called beforehand.
+ PetscScalar operator[] (PetscInt idx) const
{
- double value;
- VecGetValues(_v, 1, &idx, &value);
- return value;
+ return get(idx);
}
/*!
Get global vector
\param u Array to store the global vector. Memory allocation is needed in advance
*/
- void getGlobalVector(PetscScalar u[]);
+ void getGlobalVector(std::vector<PetscScalar>& u) const;
- /*!
- Copy local entries including ghost ones to an array
- \param u Preallocated vector for the values of local entries.
+ /* Get an entry value. This is an expensive operation,
+ and it only get local value. Use it for only test purpose
+ Get the value of an entry by [] operator.
+ setLocalAccessibleVector() must be called beforehand.
*/
- void copyValues(std::vector<double>& u) const;
-
- /// Get an entry value. This is an expensive operation,
- /// and it only get local value. Use it for only test purpose
- PetscScalar get(const PetscInt idx) const
- {
- PetscScalar x;
- VecGetValues(_v, 1, &idx, &x);
- return x;
- }
-
- // TODO eliminate in favour of getRawVector()
- /// Get PETsc vector. Use it only for test purpose
- const PETSc_Vec &getData() const
- {
- return _v;
- }
-
- // TODO preliminary
- void setZero() { VecSet(_v, 0.0); }
-
- /// Disallow moving.
- /// \todo This operator should be implemented properly when doing a
- /// general cleanup of all matrix and vector classes.
- PETScVector& operator = (PETScVector &&) = delete;
+ PetscScalar get(const PetscInt idx) const;
//! Exposes the underlying PETSc vector.
PETSc_Vec getRawVector() { return _v; }
@@ -217,8 +201,13 @@ class PETScVector
* This method is dangerous insofar as you can do arbitrary things also
* with a const PETSc vector!
*/
- PETSc_Vec getRawVector() const {return _v; }
+ PETSc_Vec getRawVector() const { return _v; }
+ /*!
+ Copy local entries including ghost ones to an array
+ \param u Preallocated vector for the values of local entries.
+ */
+ void copyValues(std::vector<PetscScalar>& u) const;
/*! View the global vector for test purpose. Do not use it for output a big vector.
\param file_name File name for output
@@ -253,7 +242,7 @@ class PETScVector
PETSc_Vec _v = nullptr;
/// Local vector, which is only for the case that _v is created
- /// with ghost entries.
+ /// with ghost entries.
mutable PETSc_Vec _v_loc = nullptr;
/// Starting index in a rank
@@ -271,19 +260,34 @@ class PETScVector
/// Flag to indicate whether the vector is created with ghost entry indices
bool _has_ghost_id = false;
+ /*!
+ \brief Array containing the entries of the vector.
+ If the vector is created without given ghost IDs, the array
+ contains all entries of the global vector, _v. Otherwise it
+ only contains the entries of the global vector owned by the
+ current rank.
+ */
+ mutable std::vector<PetscScalar> _entry_array;
+
+ /// Map global indices of ghost enrties to local indices
+ mutable std::map<PetscInt, PetscInt> _global_ids2local_ids_ghost;
+
/*!
\brief Collect local vectors
\param local_array Local array
\param global_array Global array
*/
void gatherLocalVectors(PetscScalar local_array[],
- PetscScalar global_array[]);
+ PetscScalar global_array[]) const;
/*!
Get local vector, i.e. entries in the same rank
*/
PetscScalar* getLocalVector() const;
+ /// Get local index by a global index
+ PetscInt getLocalIndex(const PetscInt global_index) const;
+
/*!
Restore array after finish access local array
\param array Pointer to the local array fetched by VecGetArray
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 19a4b01c701b95e2c59959a69346c5908c9d2a11..a365bffd7bf742bacd5883cf0e9bcdc81df08872 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -12,6 +12,7 @@
#include "MeshComponentMap.h"
+#include "BaseLib/Error.h"
#include "MeshLib/MeshSubsets.h"
#ifdef USE_PETSC
@@ -70,17 +71,18 @@ MeshComponentMap::MeshComponentMap(
for (std::size_t j = 0; j < mesh_subset.getNumberOfNodes(); j++)
{
GlobalIndexType global_id = 0;
- if (order == ComponentOrder::BY_LOCATION)
+ if (order != ComponentOrder::BY_LOCATION)
{
- global_id = static_cast<GlobalIndexType>(
- components.size() * mesh.getGlobalNodeID(j) + comp_id);
- }
- else
- {
- // _num_global_dof is used as the global index offset
- global_id = static_cast<GlobalIndexType>(
- _num_global_dof + mesh.getGlobalNodeID(j));
+ // Deactivated since this case is not suitable to
+ // arrange non ghost entries of a partition within
+ // a rank in the parallel computing.
+ OGS_FATAL("Global index in the system of equations"
+ " can only be numbered by the oder type"
+ " of ComponentOrder::BY_LOCATION");
}
+ global_id = static_cast<GlobalIndexType>(
+ components.size() * mesh.getGlobalNodeID(j)
+ + comp_id);
const bool is_ghost =
mesh.isGhostNode(mesh.getNode(j)->getID());
if (is_ghost)
@@ -294,16 +296,19 @@ GlobalIndexType MeshComponentMap::getLocalIndex(
// A special case for a ghost location with global index equal to the size
// of the local vector:
- if (-global_index == static_cast<GlobalIndexType>(_num_global_dof))
- return 0;
+ GlobalIndexType const real_global_index =
+ (-global_index == static_cast<GlobalIndexType>(_num_global_dof))
+ ? 0 : -global_index;
// TODO Find in ghost indices is O(n^2/2) for n being the length of
// _ghosts_indices. Providing an inverted table would be faster.
auto const ghost_index_it = std::find(_ghosts_indices.begin(),
- _ghosts_indices.end(), -global_index);
+ _ghosts_indices.end(),
+ real_global_index);
if (ghost_index_it == _ghosts_indices.end())
{
- OGS_FATAL("index %d not found in ghost_indices", -global_index);
+ OGS_FATAL("index %d not found in ghost_indices",
+ real_global_index);
}
// Using std::distance on a std::vector is O(1). As long as _ghost_indices
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
index 9bc318cd9d80de09761722d3f04bb93ee07de42b..5c5a15784ef8366cce321e621d41ca8166a157e9 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
@@ -27,13 +27,14 @@ LocalLinearLeastSquaresExtrapolator::LocalLinearLeastSquaresExtrapolator(
MathLib::MatrixSpecifications(dof_table.dofSizeWithoutGhosts(),
dof_table.dofSizeWithoutGhosts(),
&dof_table.getGhostIndices(),
- nullptr)))
+ nullptr))),
#ifndef USE_PETSC
- , _residuals(dof_table.size())
+ _residuals(dof_table.size()),
#else
- , _residuals(dof_table.size(), false)
+ _residuals(dof_table.dofSizeWithoutGhosts(),
+ dof_table.getGhostIndices(), false),
#endif
- , _local_to_global(dof_table)
+ _local_to_global(dof_table)
{
/* Note in case the following assertion fails:
* If you copied the extrapolation code, for your processes from
diff --git a/ProcessLib/GlobalVectorFromNamedFunction.cpp b/ProcessLib/GlobalVectorFromNamedFunction.cpp
index 6affaad855d10abcdbc36ffa7739af5d21408cfb..82020a384e3361a22acf8eed09bd3cee2ee338d2 100644
--- a/ProcessLib/GlobalVectorFromNamedFunction.cpp
+++ b/ProcessLib/GlobalVectorFromNamedFunction.cpp
@@ -8,6 +8,7 @@
*/
#include "GlobalVectorFromNamedFunction.h"
+#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/DOF/DOFTableUtil.h"
@@ -53,10 +54,10 @@ GlobalVector const& GlobalVectorFromNamedFunction::call(
_context.index = node_id;
auto const value = _function_caller.call(args);
- // TODO Problems with PETSc? (local vs. global index)
result->set(node_id, value);
}
+ MathLib::finalizeVectorAssembly(*result);
return *result;
}
diff --git a/ProcessLib/HT/Tests.cmake b/ProcessLib/HT/Tests.cmake
index 2618492602de3f525315dc2cab07980c11289eaa..c84d4ad8eba088bb98585893b92a296448f8fa08 100644
--- a/ProcessLib/HT/Tests.cmake
+++ b/ProcessLib/HT/Tests.cmake
@@ -11,3 +11,24 @@ AddTest(
ThermalConvection_const_viscosity_expected.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000.000000.vtu temperature T
VIS ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000.000000.vtu
)
+
+# MPI tests
+AddTest(
+ NAME Parallel_LARGE_2D_ThermalConvection_constviscosity
+ PATH Parabolic/HT/ConstViscosity
+ EXECUTABLE_ARGS square_5500x5500.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 4
+ TESTER vtkdiff
+ REQUIREMENTS OGS_USE_MPI
+ ABSTOL 1e-15 RELTOL 1e-14
+ DIFF_DATA
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu p p
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu p p
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu p p
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu p p
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu T T
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu T T
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu T T
+ ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu T T
+)
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 5bf2b92329e7f21b0805d592bdfd32a9b434268d..25e265394bbc0180413bf33b88b40bec6bbc273e 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -122,6 +122,8 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
StaggeredCouplingTerm const& coupling_term)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
assembleConcreteProcess(t, x, M, K, b, coupling_term);
_boundary_conditions.applyNaturalBC(t, x, K, b);
@@ -134,6 +136,9 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
GlobalVector& b, GlobalMatrix& Jac,
StaggeredCouplingTerm const& coupling_term)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+ MathLib::LinAlg::setLocalAccessibleVector(xdot);
+
assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b,
Jac, coupling_term);
@@ -235,12 +240,14 @@ void Process::computeSparsityPattern()
void Process::preTimestep(GlobalVector const& x, const double t,
const double delta_t)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
preTimestepConcreteProcess(x, t, delta_t);
_boundary_conditions.preTimestep(t);
}
void Process::postTimestep(GlobalVector const& x)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
postTimestepConcreteProcess(x);
}
@@ -248,6 +255,8 @@ void Process::computeSecondaryVariable(const double t, GlobalVector const& x,
StaggeredCouplingTerm const&
coupled_term)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
computeSecondaryVariableConcrete(t, x, coupled_term);
}
@@ -258,11 +267,13 @@ void Process::preIteration(const unsigned iter, const GlobalVector &x)
cached_var->expire();
}
+ MathLib::LinAlg::setLocalAccessibleVector(x);
preIterationConcreteProcess(iter, x);
}
NumLib::IterationResult Process::postIteration(const GlobalVector &x)
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
return postIterationConcreteProcess(x);
}
diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp
index a6c7b0da65dcde82471730a9aaa7159c9c0a0fa6..7047b87f8298689a8a914c61379bfca67d243cfe 100644
--- a/ProcessLib/StaggeredCouplingTerm.cpp
+++ b/ProcessLib/StaggeredCouplingTerm.cpp
@@ -11,10 +11,40 @@
*/
#include "StaggeredCouplingTerm.h"
+
+#include "MathLib/LinAlg/LinAlg.h"
#include "Process.h"
namespace ProcessLib
{
+
+StaggeredCouplingTerm::StaggeredCouplingTerm(
+ std::unordered_map<std::type_index, Process const&> const&
+ coupled_processes_,
+ std::unordered_map<std::type_index, GlobalVector const&> const& coupled_xs_,
+ const double dt_, const bool empty_)
+ : coupled_processes(coupled_processes_),
+ coupled_xs(coupled_xs_),
+ dt(dt_),
+ empty(empty_)
+{
+ for (auto const& coupled_x_pair : coupled_xs)
+ {
+ auto const& coupled_x = coupled_x_pair.second;
+ MathLib::LinAlg::setLocalAccessibleVector(coupled_x);
+ }
+
+ for (auto const& coupled_process_pair : coupled_processes)
+ {
+ auto const& coupled_pcs = coupled_process_pair.second;
+ auto const prevous_time_x = coupled_pcs.getPreviousTimeStepSolution();
+ if (prevous_time_x)
+ {
+ MathLib::LinAlg::setLocalAccessibleVector(*prevous_time_x);
+ }
+ }
+}
+
const StaggeredCouplingTerm createVoidStaggeredCouplingTerm()
{
std::unordered_map<std::type_index, Process const&> coupled_processes;
diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h
index 8e362536a174a17e2f3df526815b912eda275757..6cbebcc99ab36eb0feaafacb63fdd8eca9df02d6 100644
--- a/ProcessLib/StaggeredCouplingTerm.h
+++ b/ProcessLib/StaggeredCouplingTerm.h
@@ -36,13 +36,7 @@ struct StaggeredCouplingTerm
coupled_processes_,
std::unordered_map<std::type_index, GlobalVector const&> const&
coupled_xs_,
- const double dt_, const bool empty_ = false)
- : coupled_processes(coupled_processes_),
- coupled_xs(coupled_xs_),
- dt(dt_),
- empty(empty_)
- {
- }
+ const double dt_, const bool empty_ = false);
/// References to the coupled processes are distinguished by the keys of
/// process types.
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 55d9b8566ac035c443a770428c6164ada762ebae..db778b975b82bc62ee2f695cc836e2a551a9153e 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -285,6 +285,8 @@ NumLib::IterationResult TESProcess::postIterationConcreteProcess(
std::vector<double> local_x_cache;
std::vector<double> local_x_prev_ts_cache;
+ MathLib::LinAlg::setLocalAccessibleVector(*_x_previous_timestep);
+
auto check_variable_bounds = [&](std::size_t id,
TESLocalAssemblerInterface& loc_asm) {
auto const r_c_indices = NumLib::getRowColumnIndices(
@@ -340,7 +342,6 @@ TESProcess::computeVapourPartialPressure(
auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
x_mV, _assembly_params.M_react, _assembly_params.M_inert);
- // TODO Problems with PETSc? (local vs. global index)
result_cache->set(node_id, p * x_nV);
}
@@ -378,7 +379,6 @@ TESProcess::computeRelativeHumidity(
auto const p_S =
Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(T);
- // TODO Problems with PETSc? (local vs. global index)
result_cache->set(node_id, p * x_nV / p_S);
}
@@ -419,7 +419,6 @@ TESProcess::computeEquilibriumLoading(
: _assembly_params.react_sys->getEquilibriumLoading(
p_V, T, _assembly_params.M_react);
- // TODO Problems with PETSc? (local vs. global index)
result_cache->set(node_id, C_eq);
}
diff --git a/Tests/Data b/Tests/Data
index adf88547b78bf4a38568c1efbab5e848ca2911e3..92e79f16e56e4051d0d3844df9854466c2cea921 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit adf88547b78bf4a38568c1efbab5e848ca2911e3
+Subproject commit 92e79f16e56e4051d0d3844df9854466c2cea921
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 360ceed95ed88ec840d8f5a04936c74fcc9f9b4f..51582e59f1ad9f20857c41479910b892d751978f 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -95,6 +95,7 @@ void checkGlobalVectorInterfacePETSc()
const int r0 = x.getRangeBegin();
//x.get(0) is expensive, only get local value. Use it for test purpose
+ x.setLocalAccessibleVector();
ASSERT_EQ(.0, x.get(r0));
set(x, 10.);
@@ -103,18 +104,22 @@ void checkGlobalVectorInterfacePETSc()
// Value of x is not copied to y
const bool deep_copy = false;
T_VECTOR y(x, deep_copy);
+ y.setLocalAccessibleVector();
ASSERT_EQ(0, y.get(r0));
set(y, 10.0);
+ y.setLocalAccessibleVector();
ASSERT_EQ(10, y.get(r0));
// y += x
axpy(y, 1., x);
+ y.setLocalAccessibleVector();
ASSERT_EQ(20, y.get(r0));
ASSERT_EQ(80., norm2(y));
// y -= x
axpy(y, -1., x);
+ y.setLocalAccessibleVector();
ASSERT_EQ(10, y.get(r0));
ASSERT_EQ(40., norm2(y));
@@ -130,7 +135,7 @@ void checkGlobalVectorInterfacePETSc()
ASSERT_NEAR(0.0, normy-norm2(y), 1.e-10);
- double x0[16];
+ std::vector<double> x0(16);
double z[] =
{
2.0000000000000000e+01,
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 358cea04d444e8a6f576a840da55e1f9fe8cb46c..f1f1e92afa2923d5a37f734e18f69c7173b7e82c 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -166,8 +166,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
local_vec[1] = 2. * (mrank+1);
x.set(row_pos, local_vec);
- double x0[6];
- double x1[6];
+ std::vector<double> x0(6);
x.getGlobalVector(x0);
MathLib::LinAlg::matMult(A, x, b);
@@ -194,6 +193,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b,
EXPECT_GT(ls.getNumberOfIterations(), 0u);
+ std::vector<double> x1(6);
x.getGlobalVector(x1);
ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);
}
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 0350116e05203811c228d3dba7cec79450ba5c6e..74515c9b33b49e009e54e9a916d9c1878b38ef03 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -113,6 +113,7 @@ public:
) override
{
MathLib::setMatrix(M, {1.0});
+ MathLib::LinAlg::setLocalAccessibleVector(x);
MathLib::setMatrix(K, {x[0]});
MathLib::setVector(b, {0.0});
}
@@ -200,6 +201,7 @@ public:
ProcessLib::StaggeredCouplingTerm const& /*coupling_term*/
) override
{
+ MathLib::LinAlg::setLocalAccessibleVector(x_curr);
auto const u = x_curr[0];
auto const v = x_curr[1];
@@ -217,6 +219,8 @@ public:
ProcessLib::StaggeredCouplingTerm const&
coupling_term) override
{
+ MathLib::LinAlg::setLocalAccessibleVector(x_curr);
+ MathLib::LinAlg::setLocalAccessibleVector(xdot);
assemble(t, x_curr, M, K, b, coupling_term);
auto const u = x_curr[0];
diff --git a/Tests/NumLib/TestODEInt.cpp b/Tests/NumLib/TestODEInt.cpp
index 9e8ee9b0ce6b686f35e359c79893f4fe2f5da141..229f5fe7be7c328e70c5f64b586f40fd36020046 100644
--- a/Tests/NumLib/TestODEInt.cpp
+++ b/Tests/NumLib/TestODEInt.cpp
@@ -19,6 +19,7 @@
#include "BaseLib/BuildInfo.h"
#include "NumLib/NumericsConfig.h"
#include "NumLib/ODESolver/ConvergenceCriterionDeltaX.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "Tests/TestTools.h"
#include "ODEs.h"
@@ -103,6 +104,9 @@ public:
if (num_timesteps > 0)
EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb));
+ for (auto& x : sol.solutions)
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
return sol;
}
@@ -258,7 +262,8 @@ public:
TestParams::tol_picard_newton);
auto const t = sol_picard.ts[i];
- auto const sol_analyt = ODETraits<ODE>::solution(t);
+ auto sol_analyt = ODETraits<ODE>::solution(t);
+ MathLib::LinAlg::setLocalAccessibleVector(sol_analyt);
EXPECT_NEAR(sol_picard.solutions[i][comp], sol_analyt[comp],
TestParams::tol_analyt);
diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h
index 531ef5432da06a84d30d808f993998ead616d8e2..1458f420a0166eb0240ce02a0b51d50ff5ab0ff2 100644
--- a/Tests/NumLib/TimeLoopSingleODE.h
+++ b/Tests/NumLib/TimeLoopSingleODE.h
@@ -12,6 +12,7 @@
#include "NumLib/DOF/GlobalMatrixProviders.h"
#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
#include "NumLib/ODESolver/NonlinearSolver.h"
+#include "MathLib/LinAlg/LinAlg.h"
#include "ProcessLib/StaggeredCouplingTerm.h"
@@ -97,6 +98,7 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
if (time_disc.needsPreload())
{
+ MathLib::LinAlg::setLocalAccessibleVector(x);
_nonlinear_solver->assemble(x,
ProcessLib::createVoidStaggeredCouplingTerm());
time_disc.pushState(t0, x0,