From 67ed14097fd9d3ca5c2ddc839b0448ad7e718ae8 Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Wed, 1 Apr 2020 12:41:24 +0200
Subject: [PATCH] Updated Value_nz for all DCPSE operators
Signed-off-by: absingh <absingh@mpi-cbg.de>
---
src/DCPSE_op/DCPSE_Solver.hpp | 21 ++++++------
src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp | 4 +--
src/DCPSE_op/DCPSE_op.hpp | 40 ++++++++++++++++++-----
src/DCPSE_op/DCPSE_op_test2.cpp | 4 +--
4 files changed, 47 insertions(+), 22 deletions(-)
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 11618181..fad1e0a4 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -224,12 +224,12 @@ class DCPSE_scheme: public MatrixAssembler
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the term B and the Matrix A for Ax=B must contain the same number of rows\n";
return;
}
-
- if (row_b != p_map.size_local() * Sys_eqs::nvar)
- {
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set " << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar << std::endl;
- return;
- }
+ if (row_b != p_map.size_local() * Sys_eqs::nvar) {
+ std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set "
+ << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar
+ << std::endl;
+ return;
+ }
// Indicate all the non zero rows
openfpm::vector<unsigned char> nz_rows;
@@ -443,7 +443,6 @@ public:
A.resize(tot * Sys_eqs::nvar, tot * Sys_eqs::nvar,
p_map.size_local() * Sys_eqs::nvar,
p_map.size_local() * Sys_eqs::nvar);
-// std::cout << Eigen::MatrixXd(A) << std::endl;
} else
{
auto & v_cl = create_vcluster();
@@ -483,6 +482,7 @@ public:
}
+
return A;
}
@@ -494,9 +494,9 @@ public:
*/
typename Sys_eqs::Vector_type & getB(options_solver opt = options_solver::STANDARD)
{
-#ifdef SE_CLASS1
- consistency();
-#endif
+//#ifdef SE_CLASS1
+// consistency();
+//#endif
if (opt == options_solver::LAGRANGE_MULTIPLIER)
{
auto & v_cl = create_vcluster();
@@ -506,6 +506,7 @@ public:
}
}
+
return b;
}
diff --git a/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp b/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
index c9390a9d..8f6cf6c6 100644
--- a/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
+++ b/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
@@ -124,9 +124,9 @@ void MonomialBasis<dim>::generateBasis(std::vector<unsigned int> m, unsigned int
grid_key_dx_iterator_sub_bc<dim> it(grid, start, stop, bc);
// Finally compute alpha_min
- unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
+ //unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
//std::cout<<"AlphaMin: "<<alphaMin<<std::endl;
- //unsigned char alphaMin = 0; // we want to always have 1 in the basis
+ unsigned char alphaMin = 0; // we want to always have 1 in the basis
while (it.isNext())
{
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index f61c9dea..c538ffb4 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -65,9 +65,15 @@ public:
{
auto coeff_dc = dcp.getCoeffNN(key,j);
auto k = dcp.getIndexNN(key,j);
- cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
- cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp.getSign() * coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
+ auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+ auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+ //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
+
+ //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp.getSign() * coeff_dc * coeff / dcp.getEpsilonPrefactor(key);
}
}
@@ -156,10 +162,16 @@ public:
auto coeff_dc = dcp[i].getCoeffNN(key,j);
auto k = dcp[i].getIndexNN(key,j);
+ auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
- cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
- cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+ //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+ //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
}
}
}
@@ -249,10 +261,16 @@ public:
auto coeff_dc = dcp[i].getCoeffNN(key,j);
auto k = dcp[i].getIndexNN(key,j);
+ auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+ auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
- cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
- cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+ //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
}
}
}
@@ -490,10 +508,16 @@ public:
auto coeff_dc = dcp[i].getCoeffNN(key,j);
auto k = dcp[i].getIndexNN(key,j);
+ auto k_coeff = coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
- cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ auto kk_coeff = dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ o1.template value_nz<Sys_eqs>(p_map,k,cols,k_coeff,comp);
+
+
+ //cols[p_map. template getProp<0>(k)*Sys_eqs::nvar + comp] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
- cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+ //cols[p_map. template getProp<0>(key)*Sys_eqs::nvar + comp] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
}
}
}
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index dd5e6c0a..2d5f101c 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -188,10 +188,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Divergence Div(Particles, 2, rCut, 1.9);
int n=10;
double nu=1e-2;
- for(int i=0; i<=n ;i++)
+ Particles.write_frame("Stokes",0);
+ for(int i=1; i<=n ;i++)
{ RHS=-Grad(P);
DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles);
-// auto Stokes = Adv(V,V_star)-nu*Lap(V_star);
auto Stokes1 = Adv(V[0],V_star[0])-nu*Lap(V_star[0]);
auto Stokes2 = Adv(V[1],V_star[1])-nu*Lap(V_star[1]);
Solver.impose(Stokes1,bulk,RHS[0],vx);
--
GitLab