From 259b28381b6c85e16b9e14e4457940b6478ae15c Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Fri, 24 Apr 2020 17:09:10 +0200
Subject: [PATCH] Fixing DCPSE_Solver on decoupled system of equation
---
mpi_include | 2 +-
mpi_libs | 2 +-
src/DCPSE_op/DCPSEActiveGelTest.cpp | 4 +-
src/DCPSE_op/DCPSE_Solver.hpp | 56 +++++++++++++------
.../Vector/vector_dist_operators.hpp | 21 +++++--
5 files changed, 58 insertions(+), 27 deletions(-)
diff --git a/mpi_include b/mpi_include
index 16981421..4439b1ca 100644
--- a/mpi_include
+++ b/mpi_include
@@ -1 +1 @@
--I
\ No newline at end of file
+-I/home/i-bird/MPI/include
\ No newline at end of file
diff --git a/mpi_libs b/mpi_libs
index 0519ecba..3a37a5cd 100644
--- a/mpi_libs
+++ b/mpi_libs
@@ -1 +1 @@
-
\ No newline at end of file
+-Wl,-rpath -Wl,/home/i-bird/MPI/lib -Wl,--enable-new-dtags -pthread /home/i-bird/MPI/lib/libmpi.so
\ No newline at end of file
diff --git a/src/DCPSE_op/DCPSEActiveGelTest.cpp b/src/DCPSE_op/DCPSEActiveGelTest.cpp
index 7a021068..f315b46f 100644
--- a/src/DCPSE_op/DCPSEActiveGelTest.cpp
+++ b/src/DCPSE_op/DCPSEActiveGelTest.cpp
@@ -556,7 +556,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y]));
Particles.ghost_get<MolField>();
-
f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
(Pol[x] * Pol[x] + Pol[y] * Pol[y]);
@@ -1809,9 +1808,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto Pol = getV<Polarization>(Particles);
auto V = getV<Velocity>(Particles);
+ V.setVarId(0);
auto W = getV<Vorticity>(Particles);
auto g = getV<ExtForce>(Particles);
auto P = getV<Pressure>(Particles);
+ P.setVarId(2);
auto u = getV<Strain_rate>(Particles);
auto sigma = getV<Stress>(Particles);
auto h = getV<MolField>(Particles);
@@ -2041,6 +2042,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
RHS[y] = dV[y];
Particles.ghost_get<10>();
DCPSE_scheme<equations2d3pE, decltype(Particles)> Solver(Particles);
+ //Solver.setEquationStructure({eq_struct::VECTOR,eq_struct::SCALAR});
Solver.impose(Stokes1, bulk, RHS[0], vx);
Solver.impose(Stokes2, bulk, RHS[1], vy);
Solver.impose(continuity, bulkF, 0, ic);
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 0b5a4545..2bbaade2 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -43,6 +43,12 @@ enum options_solver
LAGRANGE_MULTIPLIER
};
+/*enum eq_struct
+{
+ VECTOR,
+ SCALAR
+};*/
+
//template<unsigned int prp_id> using prop_id = boost::mpl::int_<prp_id>;
template<typename Sys_eqs, typename particles_type>
@@ -70,6 +76,9 @@ class DCPSE_scheme: public MatrixAssembler
//! Particles used to impose the system
particles_type & parts;
+ //! colums shift map
+ //int col_sm[Sys_eqs::nvar];
+
//! Each point in the grid has a global id, to decompose correctly the Matrix each processor contain a
//! contiguos range of global id, example processor 0 can have from 0 to 234 and processor 1 from 235 to 512
//! no processors can have holes in the sequence, this number indicate where the sequence start for this
@@ -329,31 +338,36 @@ class DCPSE_scheme: public MatrixAssembler
public:
- /*! \brief Solve an equation
+ /*! \brief Set the structure of the system of equation
*
- * \warning exp must be a scalar type
+ * For example for stokes-flow where you are solving for V = velocity (Vector) and P = pressure (scalar)
*
- * \param exp where to store the result
+ * you should call this function with
+ *
+ * setEquationStructure({eq_struct::VECTOR,eq_struct::SCALAR})
*
*/
-/* template<typename expr_type>
- void solve(expr_type exp)
+/* void setEquationStructure(std::initializer_list<eq_struct> l)
{
- umfpack_solver<double> solver;
- auto x = solver.solve(getA(opt),getB(opt));
-
- auto parts = exp.getVector();
-
- auto it = parts.getDomainIterator();
-
- while (it.isNext())
- {
- auto p = it.get();
- exp.value(p) = x(p.getKey());
- ++it;
- }
+ int i = 0;
+ for (eq_struct e : l)
+ {
+ if (e == eq_struct::VECTOR)
+ {
+ for (int j = 0 ; j < Sys_eqs::dims ; j++)
+ {
+ col_sm[i+j] = i;
+ }
+ i += Sys_eqs::dims;
+ }
+ else
+ {
+ col_sm[i] = i;
+ }
+ }
}*/
+
/*! \brief Solve an equation
*
* \warning exp must be a scalar type
@@ -681,6 +695,12 @@ public:
if (trpl.last().row() == trpl.last().col())
is_diag = true;
+ if (trpl.last().col() == 1323)
+ {
+ int debug = 0;
+ debug++;
+ }
+
// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
}
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index b981d6fe..203656ea 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -658,6 +658,13 @@ public:
//! Property id of the point
static const unsigned int prop = prp;
+ int var_id = 0;
+
+ void setVarId(int var_id)
+ {
+ this->var_id = var_id;
+ }
+
//! constructor for an external vector
vector_dist_expression(vector & v)
:v(v)
@@ -783,7 +790,7 @@ public:
template<typename Sys_eqs, typename pmap_type, typename unordered_map_type, typename coeff_type>
inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff, unsigned int comp) const
{
- cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += coeff;
+ cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + var_id + comp] += coeff;
}
inline vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> operator[](int comp)
@@ -792,7 +799,7 @@ public:
comp_n[0] = comp;
- vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> v_exp(*this,comp_n);
+ vector_dist_expression_op<vector_dist_expression<prp,vector>,boost::mpl::int_<1>,VECT_COMP> v_exp(*this,comp_n,var_id);
return v_exp;
}
@@ -914,6 +921,8 @@ class vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP>
//! component
int comp[n];
+ int var_id = 0;
+
typedef vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP> myself;
public:
@@ -922,8 +931,8 @@ public:
//! constructor from an expresssion
- vector_dist_expression_op(const exp1 & o1, int (& comp)[n])
- :o1(o1)
+ vector_dist_expression_op(const exp1 & o1, int (& comp)[n], int var_id)
+ :o1(o1),var_id(var_id)
{
for (int i = 0 ; i < n ; i++)
{this->comp[i] = comp[i];}
@@ -1004,7 +1013,7 @@ public:
#endif
- o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp_ + comp[0]);
+ o1.template value_nz<Sys_eqs>(p_map,key,cols,coeff,comp_ + var_id + comp[0]);
}
inline vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> operator[](int comp_)
@@ -1015,7 +1024,7 @@ public:
{comp_n[i] = comp[i];}
comp_n[n] = comp_;
- vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> v_exp(o1,comp_n);
+ vector_dist_expression_op<exp1,boost::mpl::int_<2>,VECT_COMP> v_exp(o1,comp_n,var_id);
return v_exp;
}
--
GitLab