COOLFluiD  Release kernel
COOLFluiD is a Collaborative Simulation Environment (CSE) focused on complex MultiPhysics simulations.
utest-lss-atomic.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010-2013 von Karman Institute for Fluid Dynamics, Belgium
2 //
3 // This software is distributed under the terms of the
4 // GNU Lesser General Public License version 3 (LGPLv3).
5 // See doc/lgpl.txt and doc/gpl.txt for the license text.
6 //
7 
8 #define BOOST_TEST_DYN_LINK
9 #define BOOST_TEST_MODULE "Test module for cf3::math::LSS where testing indivdual operations."
10 
12 
13 #include <fstream>
14 
15 #include <boost/test/unit_test.hpp>
16 #include <boost/assign/std/vector.hpp>
17 #include <boost/lexical_cast.hpp>
18 
19 #include "common/Log.hpp"
20 #include "math/LSS/System.hpp"
21 #include "math/VariablesDescriptor.hpp"
22 
24 #include "common/PE/debug.hpp"
25 #include "common/Environment.hpp"
26 
28 
29 using namespace boost::assign;
30 
31 using namespace cf3;
32 using namespace cf3::math;
33 using namespace cf3::math::LSS;
34 
36 
37 struct LSSAtomicFixture
38 {
40  LSSAtomicFixture() :
41  solvertype("Trilinos"),
42  gid(0),
43  irank(0),
44  nproc(1),
45  neq(2),
46  rank_updatable(0),
47  node_connectivity(0),
48  starting_indices(0),
49  blockcol_size(0),
50  blockrow_size(0)
51  {
52  if (common::PE::Comm::instance().is_initialized())
53  {
54  nproc=common::PE::Comm::instance().size();
55  irank=common::PE::Comm::instance().rank();
56  BOOST_CHECK_EQUAL(nproc,2);
57  }
58  m_argc = boost::unit_test::framework::master_test_suite().argc;
59  m_argv = boost::unit_test::framework::master_test_suite().argv;
60 
61  if(m_argc != 2)
62  throw common::ParsingFailed(FromHere(), "Failed to parse command line arguments: expected one argument: builder name for the matrix");
63  matrix_builder = m_argv[1];
64  }
65 
67  ~LSSAtomicFixture()
68  {
69  }
70 
72  void build_commpattern(common::PE::CommPattern& cp)
73  {
74  if (irank==0)
75  {
76  gid += 1,2,8,7,3,4,5,6;
77  rank_updatable += 0,1,0,0,1,0,0,1;
78  } else {
79  gid += 5,0,2,7,1,3,8,4,6;
80  rank_updatable += 0,1,1,0,0,1,0,0,1;
81  }
82  cp.insert("gid",gid,1,false);
83  cp.setup(Handle<common::PE::CommWrapper>(cp.get_child("gid")),rank_updatable);
84  }
85 
87  void build_system(LSS::System& sys, common::PE::CommPattern& cp)
88  {
89  starting_indices.clear();
90  node_connectivity.clear();
91  if (irank==0)
92  {
93  node_connectivity += 0,2,4,6,1,2,3,5,1,3,5,7,1,2,3,5,0,1,3,4,5,6;
94  starting_indices += 0,4,4,8,12,12,16,22,22;
95  blockcol_size = 8;
96  blockrow_size = 5;
97  } else {
98  node_connectivity += 1,2,5,2,5,8,0,7,3,1,2,4,5,6,7,8,2,5,8;
99  starting_indices += 0,0,3,9,9,9,16,16,16,19;
100  blockcol_size = 9;
101  blockrow_size = 4;
102  }
103  sys.create(cp,neq,node_connectivity,starting_indices);
104  }
105 
107  std::string solvertype;
108  std::string matrix_builder;
109 
111  int irank;
112  int nproc;
113  int neq;
114  int m_argc;
115  char** m_argv;
116 
118  std::vector<Uint> gid;
119  std::vector<Uint> rank_updatable;
120 
122  std::vector<Uint> node_connectivity;
123  std::vector<Uint> starting_indices;
124  int blockrow_size;
125  int blockcol_size;
126 
127 };
128 
130 
131 BOOST_FIXTURE_TEST_SUITE( LSSAtomicSuite, LSSAtomicFixture )
132 
133 
135 BOOST_AUTO_TEST_CASE( init_mpi )
136 {
137  common::PE::Comm::instance().init(m_argc,m_argv);
138  BOOST_CHECK_EQUAL(common::PE::Comm::instance().is_active(),true);
139  CFinfo.setFilterRankZero(false);
140  common::Core::instance().environment().options().set("log_level", 4u);
141  common::Core::instance().environment().options().set("exception_backtrace", false);
142  common::Core::instance().environment().options().set("exception_outputs", false);
143 // common::PE::wait_for_debugger(1);
144 }
145 
147 
148 BOOST_AUTO_TEST_CASE( test_matrix_only )
149 {
150 
151  // build a commpattern and a matrix
152  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
153  common::PE::CommPattern& cp = *cp_ptr;
154  build_commpattern(cp);
155  boost::shared_ptr<LSS::System> sys(common::allocate_component<LSS::System>("sys"));
156  sys->options().option("matrix_builder").change_value(matrix_builder);
157  build_system(*sys,cp);
158  Handle<LSS::Matrix> mat=sys->matrix();
159  BOOST_CHECK_EQUAL(mat->is_created(),true);
160  BOOST_CHECK_EQUAL(mat->solvertype(),solvertype);
161  BOOST_CHECK_EQUAL(mat->neq(),neq);
162  BOOST_CHECK_EQUAL(mat->blockrow_size(),blockrow_size);
163  BOOST_CHECK_EQUAL(mat->blockcol_size(),blockcol_size);
164 
165  // just to see if its crashing or not
166 // mat->print(std::cout);
167 // mat->print(CFinfo);
168  mat->print("test_matrix_" + boost::lexical_cast<std::string>(irank) + ".plt");
169 
170  // counter-checking data
171  std::vector<Uint> cols(0);
172  std::vector<Uint> rows(0);
173  std::vector<Real> vals(0);
174 
175  // check reset
176  mat->reset(1.);
177  mat->debug_data(rows,cols,vals);
178  BOOST_CHECK_EQUAL(cols.size(),node_connectivity.size()*neq*neq);
179  BOOST_CHECK_EQUAL(rows.size(),node_connectivity.size()*neq*neq);
180  BOOST_CHECK_EQUAL(vals.size(),node_connectivity.size()*neq*neq);
181  BOOST_FOREACH(double v, vals) BOOST_CHECK_EQUAL(v,1.);
182  mat->reset(0.);
183  mat->debug_data(rows,cols,vals);
184  BOOST_FOREACH(double v, vals) BOOST_CHECK_EQUAL(v,0.);
185 
186  // diagonal access check
187  mat->reset();
188  std::vector<Real> diag(blockcol_size*neq);
189  for (int i=0; i<diag.size(); i++) diag[i]=i;
190  mat->set_diagonal(diag);
191  mat->add_diagonal(diag);
192  diag.clear();
193  mat->get_diagonal(diag);
194  BOOST_CHECK_EQUAL(diag.size(),blockcol_size*neq);
195  mat->debug_data(rows,cols,vals);
196  for (int i=0; i<(const int)vals.size(); i++)
197  {
198  if (cp.isUpdatable()[rows[i]/neq]) { BOOST_CHECK_EQUAL(diag[rows[i]],2.*rows[i]); }
199  else { BOOST_CHECK_EQUAL(diag[rows[i]],0.); }
200  if (rows[i]==cols[i]) { BOOST_CHECK_EQUAL(vals[i],diag[rows[i]]); }
201  else { BOOST_CHECK_EQUAL(vals[i],0.); }
202  }
203 
204  // individual access
205  mat->reset();
206  if (irank==0)
207  {
208  mat->set_value(6,6,1.);
209  mat->set_value(7,6,1.);
210  mat->add_value(7,6,1.);
211  mat->set_value(6,7,3.);
212  mat->add_value(7,7,4.);
213  Real v;
214  mat->get_value(6,6,v);
215  BOOST_CHECK_EQUAL(v,1.);
216  mat->get_value(7,6,v);
217  BOOST_CHECK_EQUAL(v,2.);
218  mat->get_value(6,7,v);
219  BOOST_CHECK_EQUAL(v,3.);
220  mat->get_value(7,7,v);
221  BOOST_CHECK_EQUAL(v,4.);
222  mat->debug_data(rows,cols,vals);
223  for (int i=0; i<(const int)vals.size(); i++)
224  if ((rows[i]<6)||(rows[i]>7))
225  BOOST_CHECK_EQUAL(vals[i],0.);
226  for (int i=0; i<(const int)vals.size(); i++)
227  if ((cols[i]<6)||(cols[i]>7))
228  BOOST_CHECK_EQUAL(vals[i],0.);
229  }
230 
231  // performant access
232  mat->reset();
233  if (irank==1)
234  {
235  LSS::BlockAccumulator ba;
236  ba.resize(3,neq);
237  ba.mat << 53., 54., 51., 52., 55., 56.,
238  59., 60., 57., 58., 61., 62.,
239  23., 24., 21., 22., 25., 26.,
240  29., 30., 27., 28., 31., 32.,
241  83., 84., 81., 82., 85., 86.,
242  89., 90., 87., 88., 91., 92.;
243  ba.indices[0]=5;
244  ba.indices[1]=2;
245  ba.indices[2]=8;
246  mat->set_values(ba);
247  mat->add_values(ba);
248  ba.reset();
249  ba.indices[0]=2;
250  ba.indices[1]=5;
251  ba.indices[2]=8;
252  mat->get_values(ba);
253 
254  mat->print("test_ba_assembly_" + boost::lexical_cast<std::string>(irank) + ".plt");
255 
256  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(0,i-1),(double)((ba.indices[0]*10+i+0)*2));
257  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(1,i-1),(double)((ba.indices[0]*10+i+6)*2));
258  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(2,i-1),(double)((ba.indices[1]*10+i+0)*2));
259  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(3,i-1),(double)((ba.indices[1]*10+i+6)*2));
260  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(4,i-1),(double)((ba.indices[2]*10+i+0)*2));
261  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(5,i-1),(double)((ba.indices[2]*10+i+6)*2));
262  mat->debug_data(rows,cols,vals);
263  Real ctr=1.;
264  for (int i=0; i<(const int)vals.size(); i++)
265  if (((rows[i]/neq==ba.indices[0])||(rows[i]/neq==ba.indices[1])||(rows[i]/neq==ba.indices[2]))&&
266  ((cols[i]/neq==ba.indices[0])||(cols[i]/neq==ba.indices[1])||(cols[i]/neq==ba.indices[2])))
267  {
268  BOOST_CHECK_EQUAL(vals[i],(double)((rows[i]/neq*10.+ctr)*2));
269  ctr+=1.;
270  if (ctr==13.) ctr=1.;
271  } else {
272  BOOST_CHECK_EQUAL(vals[i],0.);
273  }
274  }
275 
276 
277 
278  // performant access - out of range access does not fail
279  mat->reset();
280  if (irank==1)
281  {
282  LSS::BlockAccumulator ba;
283  ba.resize(3,neq);
284  ba.mat << 99., 99., 99., 99., 99., 99.,
285  99., 99., 99., 99., 99., 99.,
286  23., 24., 21., 22., 25., 26.,
287  29., 30., 27., 28., 31., 32.,
288  99., 99., 99., 99., 99., 99.,
289  99., 99., 99., 99., 99., 99.;
290  ba.indices[0]=3;
291  ba.indices[1]=2;
292  ba.indices[2]=7;
293  mat->set_values(ba);
294  mat->add_values(ba);
295  ba.reset(-1.);
296  ba.indices[0]=7;
297  ba.indices[1]=3;
298  ba.indices[2]=2;
299  mat->get_values(ba);
300  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(0,i-1),0.);
301  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(1,i-1),0.);
302  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(2,i-1),0.);
303  for (int i=1; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(3,i-1),0.);
304  for (int i=1; i<3; i++) BOOST_CHECK_EQUAL(ba.mat(4,i-1),(double)((ba.indices[2]*10+i+4+0)*2));
305  for (int i=1; i<3; i++) BOOST_CHECK_EQUAL(ba.mat(5,i-1),(double)((ba.indices[2]*10+i+4+6)*2));
306  for (int i=3; i<5; i++) BOOST_CHECK_EQUAL(ba.mat(4,i-1),(double)((ba.indices[2]*10+i+0+0)*2));
307  for (int i=3; i<5; i++) BOOST_CHECK_EQUAL(ba.mat(5,i-1),(double)((ba.indices[2]*10+i+0+6)*2));
308  for (int i=5; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(4,i-1),(double)((ba.indices[2]*10+i-4+0)*2));
309  for (int i=5; i<7; i++) BOOST_CHECK_EQUAL(ba.mat(5,i-1),(double)((ba.indices[2]*10+i-4+6)*2));
310  mat->debug_data(rows,cols,vals);
311  Real ctr=1.;
312  for (int i=0; i<(const int)vals.size(); i++)
313  if ((rows[i]/neq==ba.indices[2])&&
314  ((cols[i]/neq==ba.indices[0])||(cols[i]/neq==ba.indices[1])||(cols[i]/neq==ba.indices[2])))
315  {
316  BOOST_CHECK_EQUAL(vals[i],(double)((rows[i]/neq*10.+ctr)*2));
317  ctr+=1.;
318  if (ctr==13.) ctr=1.;
319  } else {
320  BOOST_CHECK_EQUAL(vals[i],0.);
321  }
322  }
323 
324  // bc-related: dirichlet-condition
325  mat->reset(-1.);
326  if (irank==0)
327  {
328  mat->set_row(3,1,1.,0.);
329  mat->debug_data(rows,cols,vals);
330  for (int i=0; i<(const int)vals.size(); i++)
331  {
332  if (rows[i]==7)
333  {
334  if (cols[i]==7) { BOOST_CHECK_EQUAL(vals[i],1.); }
335  else { BOOST_CHECK_EQUAL(vals[i],0.); }
336  } else {
337  BOOST_CHECK_EQUAL(vals[i],-1.);
338  }
339  }
340  }
341 
342  // bc-related: dirichlet-condition doesnt fail for ghost node
343  mat->reset(-1.);
344  if (irank==0)
345  {
346  mat->set_row(1,1,1.,0.);
347  mat->debug_data(rows,cols,vals);
348  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,-1.);
349  }
350 
351  // bc-related: symmetricizing a dirichlet
352  try
353  {
354  mat->reset(1.);
355  if (irank==0)
356  {
357  mat->set_value(11,4,5.);
358  mat->set_value(11,5,6.);
359  mat->set_value(11,6,7.);
360  mat->set_value(11,7,8.);
361  mat->set_value(11,10,11.);
362  mat->set_value(11,11,12.);
363  mat->set_value(11,12,13.);
364  mat->set_value(11,13,14.);
365  mat->get_column_and_replace_to_zero(5,1,vals);
366  vals[0]+=1.;
367  vals[1]+=2.;
368  vals[2]+=3.;
369  vals[3]+=4.;
370  vals[8]+=9.;
371  vals[9]+=10.;
372  vals[14]+=15.;
373  vals[15]+=16.;
374  for(int i=0; i<vals.size(); i++) BOOST_CHECK_EQUAL(vals[i],(double)(i+1));
375  mat->debug_data(rows,cols,vals);
376  for (int i=0; i<(const int)vals.size(); i++)
377  {
378  if (cols[i]==11)
379  {
380  BOOST_CHECK_EQUAL(vals[i],0.);
381  } else {
382  BOOST_CHECK_EQUAL(vals[i],1.);
383  }
384  }
385  }
386  }
387  catch(common::NotImplemented&)
388  {
389  CFinfo << "skipping symmetric dirichlet test" << CFendl;
390  }
391 
392  // bc-related: periodicity
393  mat->reset(-2.);
394  if (irank==0)
395  {
396  mat->set_value( 2, 4,40.);
397  mat->set_value( 3, 4,41.);
398  mat->set_value( 4, 4,21.); // half!
399  mat->set_value( 5, 4,21.5); // half!
400  mat->set_value( 6, 4,44.);
401  mat->set_value( 7, 4,45.);
402  mat->set_value(10, 4,23.); // half!
403  mat->set_value(11, 4,23.5); // half!
404  mat->set_value( 2, 5,48.);
405  mat->set_value( 3, 5,49.);
406  mat->set_value( 4, 5,25.); // half!
407  mat->set_value( 5, 5,25.5); // half!
408  mat->set_value( 6, 5,52.);
409  mat->set_value( 7, 5,53.);
410  mat->set_value(10, 5,27.); // half!
411  mat->set_value(11, 5,27.5); // half!
412 
413  mat->set_value(11,11, 5.); // half!
414  mat->set_value(10,11, 5.5); // half!
415  mat->set_value( 7,11,12.);
416  mat->set_value( 6,11,13.);
417  mat->set_value( 5,11, 7.); // half!
418  mat->set_value( 4,11, 7.5); // half!
419  mat->set_value( 3,11,16.);
420  mat->set_value( 2,11,17.);
421  mat->set_value(11,10, 9.); // half!
422  mat->set_value(10,10, 9.5); // half!
423  mat->set_value( 7,10,20.);
424  mat->set_value( 6,10,21.);
425  mat->set_value( 5,10,11.); // half!
426  mat->set_value( 4,10,11.5); // half!
427  mat->set_value( 3,10,24.);
428  mat->set_value( 2,10,25.);
429 
430  mat->tie_blockrow_pairs(2,5);
431 
432  mat->debug_data(rows,cols,vals);
433  for (int i=0; i<(const int)vals.size(); i++)
434  if ((rows[i]==4)||(rows[i]==5))
435  {
436  if ((cols[i]!=10)&&(cols[i]!=11)) { BOOST_CHECK_EQUAL(vals[i],65.); }
437  else { BOOST_CHECK_EQUAL(vals[i],0.); }
438  }
439 
440  for (int i=0; i<(const int)vals.size(); i++)
441  if (rows[i]==10)
442  {
443  if (cols[i]==4) { BOOST_CHECK_EQUAL(vals[i],-1.); }
444  else if (cols[i]==10) { BOOST_CHECK_EQUAL(vals[i],1.); }
445  else { BOOST_CHECK_EQUAL(vals[i],0.); }
446  }
447 
448  for (int i=0; i<(const int)vals.size(); i++)
449  if (rows[i]==11)
450  {
451  if (cols[i]==5) { BOOST_CHECK_EQUAL(vals[i],-1.); }
452  else if (cols[i]==11) { BOOST_CHECK_EQUAL(vals[i],1.); }
453  else { BOOST_CHECK_EQUAL(vals[i],0.); }
454  }
455 
456  for (int i=0; i<(const int)vals.size(); i++)
457  if ((rows[i]!=4)&&(rows[i]!=5)&&(rows[i]!=10)&&(rows[i]!=11))
458  BOOST_CHECK_EQUAL(vals[i],-2.);
459  }
460 
461  // bc-related: periodicity does not fail
462  mat->reset(-2.);
463  if (irank==0)
464  {
465  mat->tie_blockrow_pairs(1,4);
466  mat->debug_data(rows,cols,vals);
467  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,-2.);
468  }
469 
470  // post-destroy check
471  mat->destroy();
472  BOOST_CHECK_EQUAL(mat->is_created(),false);
473 }
474 
476 
477 BOOST_AUTO_TEST_CASE( test_vector_only )
478 {
479  // build a commpattern and the two vectors
480  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
481  common::PE::CommPattern& cp = *cp_ptr;
482  build_commpattern(cp);
483  boost::shared_ptr<LSS::System> sys(common::allocate_component<LSS::System>("sys"));
484  sys->options().option("matrix_builder").change_value(matrix_builder);
485  build_system(*sys,cp);
486  Handle<LSS::Vector> sol=sys->solution();
487  Handle<LSS::Vector> rhs=sys->rhs();
488  BOOST_CHECK_EQUAL(sol->is_created(),true);
489  BOOST_CHECK_EQUAL(rhs->is_created(),true);
490  BOOST_CHECK_EQUAL(sol->solvertype(),solvertype);
491  BOOST_CHECK_EQUAL(rhs->solvertype(),solvertype);
492  BOOST_CHECK_EQUAL(sol->neq(),neq);
493  BOOST_CHECK_EQUAL(rhs->neq(),neq);
494  BOOST_CHECK_EQUAL(sol->blockrow_size(),blockcol_size); // needs to include ghost entries
495  BOOST_CHECK_EQUAL(rhs->blockrow_size(),blockcol_size); // needs to include ghost entries
496 
497  // just to see if its crashing or not
498 // sol->print(std::cout);
499 // sol->print(CFinfo);
500  sol->print("test_vector_" + boost::lexical_cast<std::string>(irank) + ".plt");
501 
502  // counter-checking data
503  std::vector<Real> vals(0);
504  Real val;
505 
506  // check reset
507  sol->reset(1.);
508  sol->debug_data(vals);
509  BOOST_CHECK_EQUAL(vals.size(),gid.size()*neq);
510  BOOST_FOREACH(double v, vals) BOOST_CHECK_EQUAL(v,1.);
511  sol->reset();
512  sol->debug_data(vals);
513  BOOST_FOREACH(double v, vals) BOOST_CHECK_EQUAL(v,0.);
514 
515  // set by row-wise index
516  sol->reset();
517  sol->set_value(5,1.);
518  sol->add_value(5,1.);
519  val=0.;
520  sol->get_value(5,val);
521  BOOST_CHECK_EQUAL(val,2.);
522  sol->debug_data(vals);
523  for (int i=0; i<vals.size(); i++)
524  {
525  if (i==5) { BOOST_CHECK_EQUAL(vals[i],2.); }
526  else { BOOST_CHECK_EQUAL(vals[i],0.); }
527  }
528 
529  // set by block-wise index
530  sol->reset();
531  sol->set_value(2,1,1.);
532  sol->add_value(2,1,1.);
533  val=0.;
534  sol->get_value(2,1,val);
535  BOOST_CHECK_EQUAL(val,2.);
536  sol->debug_data(vals);
537  for (int i=0; i<vals.size(); i++)
538  {
539  if (i==5) { BOOST_CHECK_EQUAL(vals[i],2.); }
540  else { BOOST_CHECK_EQUAL(vals[i],0.); }
541  }
542 
543  // set via blockaccumulator
544  sol->reset();
545  BlockAccumulator ba;
546  ba.resize(3,neq);
547  ba.rhs << 8.,9. , 2.,3. , 6.,7.;
548  ba.indices[0]=4;
549  ba.indices[1]=1;
550  ba.indices[2]=3;
551  sol->set_rhs_values(ba);
552  sol->add_rhs_values(ba);
553  ba.reset();
554  sol->get_rhs_values(ba);
555  sol->debug_data(vals);
556  for (int i=0; i<vals.size(); i++)
557  {
558  if ((i/neq==ba.indices[0])||(i/neq==ba.indices[1])||(i/neq==ba.indices[2]))
559  {
560  BOOST_CHECK_EQUAL(vals[i],(double)(i*2));
561  } else {
562  BOOST_CHECK_EQUAL(vals[i],0.);
563  }
564  }
565  BOOST_CHECK_EQUAL(ba.rhs[0],16.);
566  BOOST_CHECK_EQUAL(ba.rhs[1],18.);
567  BOOST_CHECK_EQUAL(ba.rhs[2],4.);
568  BOOST_CHECK_EQUAL(ba.rhs[3],6.);
569  BOOST_CHECK_EQUAL(ba.rhs[4],12.);
570  BOOST_CHECK_EQUAL(ba.rhs[5],14.);
571  BOOST_CHECK_EQUAL(ba.sol.isConstant(0.,1.e-10),true);
572 
573  sol->reset();
574  ba.sol << 80.,90. , 20.,30. , 60.,70.;
575  sol->set_sol_values(ba);
576  sol->add_sol_values(ba);
577  ba.reset();
578  sol->get_sol_values(ba);
579  sol->debug_data(vals);
580  for (int i=0; i<vals.size(); i++)
581  {
582  if ((i/neq==ba.indices[0])||(i/neq==ba.indices[1])||(i/neq==ba.indices[2]))
583  {
584  BOOST_CHECK_EQUAL(vals[i],(double)(i*20.));
585  } else {
586  BOOST_CHECK_EQUAL(vals[i],0.);
587  }
588  }
589  BOOST_CHECK_EQUAL(ba.sol[0],160.);
590  BOOST_CHECK_EQUAL(ba.sol[1],180.);
591  BOOST_CHECK_EQUAL(ba.sol[2],40.);
592  BOOST_CHECK_EQUAL(ba.sol[3],60.);
593  BOOST_CHECK_EQUAL(ba.sol[4],120.);
594  BOOST_CHECK_EQUAL(ba.sol[5],140.);
595  BOOST_CHECK_EQUAL(ba.rhs.isConstant(0.,1.e-10),true);
596 
597  // check destroy
598  sol->destroy();
599  BOOST_CHECK_EQUAL(sol->is_created(),false);
600 }
601 
603 
604 BOOST_AUTO_TEST_CASE( test_complete_system )
605 {
606  // build a commpattern and the system
607  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
608  common::PE::CommPattern& cp = *cp_ptr;
609  build_commpattern(cp);
610  boost::shared_ptr<LSS::System> sys(common::allocate_component<LSS::System>("sys"));
611  sys->options().option("matrix_builder").change_value(matrix_builder);
612  build_system(*sys,cp);
613  BOOST_CHECK_EQUAL(sys->is_created(),true);
614  BOOST_CHECK_EQUAL(sys->solvertype(),solvertype);
615 
616  // just to see if its crashing or not
617 // sys->print(std::cout);
618 // sys->print(CFinfo);
619  sys->print("test_system_" + boost::lexical_cast<std::string>(irank) + ".plt");
620 
621  // counter-checking data
622  std::vector<Uint> cols(0);
623  std::vector<Uint> rows(0);
624  std::vector<Real> vals(0);
625 
626  // reset
627  sys->reset(1.);
628  sys->matrix()->debug_data(rows,cols,vals);
629  BOOST_CHECK_EQUAL(cols.size(),node_connectivity.size()*neq*neq);
630  BOOST_CHECK_EQUAL(rows.size(),node_connectivity.size()*neq*neq);
631  BOOST_CHECK_EQUAL(vals.size(),node_connectivity.size()*neq*neq);
632  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,1.);
633  sys->solution()->debug_data(vals);
634  BOOST_CHECK_EQUAL(vals.size(),gid.size()*neq);
635  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,1.);
636  sys->rhs()->debug_data(vals);
637  BOOST_CHECK_EQUAL(vals.size(),gid.size()*neq);
638  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,1.);
639  sys->reset();
640  sys->matrix()->debug_data(rows,cols,vals);
641  BOOST_CHECK_EQUAL(cols.size(),node_connectivity.size()*neq*neq);
642  BOOST_CHECK_EQUAL(rows.size(),node_connectivity.size()*neq*neq);
643  BOOST_CHECK_EQUAL(vals.size(),node_connectivity.size()*neq*neq);
644  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,0.);
645  sys->solution()->debug_data(vals);
646  BOOST_CHECK_EQUAL(vals.size(),gid.size()*neq);
647  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,0.);
648  sys->rhs()->debug_data(vals);
649  BOOST_CHECK_EQUAL(vals.size(),gid.size()*neq);
650  BOOST_FOREACH(Real i,vals) BOOST_CHECK_EQUAL(i,0.);
651 
652  // dirichlet bc
653  sys->matrix()->reset(2.);
654  sys->solution()->reset(3.);
655  sys->rhs()->reset(4.);
656  if (irank==0)
657  {
658  sys->dirichlet(3,1,5.,false);
659  sys->matrix()->debug_data(rows,cols,vals);
660  for (int i=0; i<vals.size(); i++)
661  {
662  if (rows[i]==7)
663  {
664  if (cols[i]==7) { BOOST_CHECK_EQUAL(vals[i],1.); }
665  else { BOOST_CHECK_EQUAL(vals[i],0.); }
666  } else {
667  if (cols[i]==7)
668  {
669  //BOOST_CHECK_EQUAL(vals[i],0.); // only works for symmetric Dirichlet BC
670  }
671  else
672  {
673  BOOST_CHECK_EQUAL(vals[i],2.);
674  }
675  }
676  }
677  sys->solution()->debug_data(vals);
678  for (int i=0; i<vals.size(); i++)
679  {
680  if (i==7) { BOOST_CHECK_EQUAL(vals[i],5.); }
681  else { BOOST_CHECK_EQUAL(vals[i],3.); }
682  }
683  sys->rhs()->debug_data(vals);
684  for (int i=0; i<4; i++) BOOST_CHECK_EQUAL(vals[i],4.);
685  for (int i=4; i<7; i++) BOOST_CHECK_EQUAL(vals[i],4.);
686  for (int i=7; i<8; i++) BOOST_CHECK_EQUAL(vals[i],5.);
687  for (int i=8; i<10; i++) BOOST_CHECK_EQUAL(vals[i],4.);
688  for (int i=10; i<14; i++) BOOST_CHECK_EQUAL(vals[i],4.);
689  for (int i=14; i<16; i++) BOOST_CHECK_EQUAL(vals[i],4.);
690  }
691 
692  // performant access - out of range access does not fail
693  sys->reset();
694  if (irank==1)
695  {
696  LSS::BlockAccumulator ba;
697  ba.resize(3,neq);
698  ba.mat << 99., 99., 99., 99., 99., 99.,
699  99., 99., 99., 99., 99., 99.,
700  1., 1., 1., 1., 1., 1.,
701  1., 1., 1., 1., 1., 1.,
702  99., 99., 99., 99., 99., 99.,
703  99., 99., 99., 99., 99., 99.;
704  ba.rhs << 2., 2., 2., 2., 2., 2.;
705  ba.sol << 3., 3., 3., 3., 3., 3.;
706  ba.indices[0]=3;
707  ba.indices[1]=2;
708  ba.indices[2]=7;
709  sys->set_values(ba);
710  sys->add_values(ba);
711  ba.reset();
712  ba.indices[0]=7;
713  ba.indices[1]=3;
714  ba.indices[2]=2;
715  sys->get_values(ba);
716  sys->matrix()->debug_data(rows,cols,vals);
717  for (int i=0; i<(const int)vals.size(); i++)
718  {
719  if ((rows[i]/neq==ba.indices[2])&&
720  ((cols[i]/neq==ba.indices[0])||(cols[i]/neq==ba.indices[1])||(cols[i]/neq==ba.indices[2]))) { BOOST_CHECK_EQUAL(vals[i],2.); }
721  else { BOOST_CHECK_EQUAL(vals[i],0.); }
722  }
723  sys->rhs()->debug_data(vals);
724  for (int i=0; i<(const int)vals.size(); i++)
725  {
726  if ((i/neq==ba.indices[0])||(i/neq==ba.indices[1])||(i/neq==ba.indices[2])) { BOOST_CHECK_EQUAL(vals[i],4.); }
727  else { BOOST_CHECK_EQUAL(vals[i],0.); }
728  }
729  sys->solution()->debug_data(vals);
730  for (int i=0; i<(const int)vals.size(); i++)
731  {
732  if ((i/neq==ba.indices[0])||(i/neq==ba.indices[1])||(i/neq==ba.indices[2])) { BOOST_CHECK_EQUAL(vals[i],6.); }
733  else { BOOST_CHECK_EQUAL(vals[i],0.); }
734  }
735  }
736 
737  // check periodicity
738  sys->matrix()->reset(-2.);
739  sys->solution()->reset(-3.);
740  sys->rhs()->reset(-4.);
741  if (irank==0)
742  {
743  sys->matrix()->set_value( 2, 4,40.);
744  sys->matrix()->set_value( 3, 4,41.);
745  sys->matrix()->set_value( 4, 4,21.); // half!
746  sys->matrix()->set_value( 5, 4,21.5); // half!
747  sys->matrix()->set_value( 6, 4,44.);
748  sys->matrix()->set_value( 7, 4,45.);
749  sys->matrix()->set_value(10, 4,23.); // half!
750  sys->matrix()->set_value(11, 4,23.5); // half!
751  sys->matrix()->set_value( 2, 5,48.);
752  sys->matrix()->set_value( 3, 5,49.);
753  sys->matrix()->set_value( 4, 5,25.); // half!
754  sys->matrix()->set_value( 5, 5,25.5); // half!
755  sys->matrix()->set_value( 6, 5,52.);
756  sys->matrix()->set_value( 7, 5,53.);
757  sys->matrix()->set_value(10, 5,27.); // half!
758  sys->matrix()->set_value(11, 5,27.5); // half!
759 
760  sys->matrix()->set_value(11,11, 5.); // half!
761  sys->matrix()->set_value(10,11, 5.5); // half!
762  sys->matrix()->set_value( 7,11,12.);
763  sys->matrix()->set_value( 6,11,13.);
764  sys->matrix()->set_value( 5,11, 7.); // half!
765  sys->matrix()->set_value( 4,11, 7.5); // half!
766  sys->matrix()->set_value( 3,11,16.);
767  sys->matrix()->set_value( 2,11,17.);
768  sys->matrix()->set_value(11,10, 9.); // half!
769  sys->matrix()->set_value(10,10, 9.5); // half!
770  sys->matrix()->set_value( 7,10,20.);
771  sys->matrix()->set_value( 6,10,21.);
772  sys->matrix()->set_value( 5,10,11.); // half!
773  sys->matrix()->set_value( 4,10,11.5); // half!
774  sys->matrix()->set_value( 3,10,24.);
775  sys->matrix()->set_value( 2,10,25.);
776 
777  sys->solution()->set_value(4,-1.);
778  sys->solution()->set_value(5, 0.);
779  sys->solution()->set_value(10,-5.);
780  sys->solution()->set_value(11,-6.);
781 
782  sys->rhs()->set_value(4,3.);
783  sys->rhs()->set_value(5,4.);
784  sys->rhs()->set_value(10,-7.);
785  sys->rhs()->set_value(11,-8.);
786 
787  sys->periodicity(2,5);
788 
789  sys->matrix()->debug_data(rows,cols,vals);
790  for (int i=0; i<(const int)vals.size(); i++)
791  if ((rows[i]==4)||(rows[i]==5))
792  {
793  if ((cols[i]!=10)&&(cols[i]!=11)) { BOOST_CHECK_EQUAL(vals[i],65.); }
794  else { BOOST_CHECK_EQUAL(vals[i],0.); }
795  }
796 
797  for (int i=0; i<(const int)vals.size(); i++)
798  if (rows[i]==10)
799  {
800  if (cols[i]==4) { BOOST_CHECK_EQUAL(vals[i],-1.); }
801  else if (cols[i]==10) { BOOST_CHECK_EQUAL(vals[i],1.); }
802  else { BOOST_CHECK_EQUAL(vals[i],0.); }
803  }
804 
805  for (int i=0; i<(const int)vals.size(); i++)
806  if (rows[i]==11)
807  {
808  if (cols[i]==5) { BOOST_CHECK_EQUAL(vals[i],-1.); }
809  else if (cols[i]==11) { BOOST_CHECK_EQUAL(vals[i],1.); }
810  else { BOOST_CHECK_EQUAL(vals[i],0.); }
811  }
812 
813  for (int i=0; i<(const int)vals.size(); i++)
814  if ((rows[i]!=4)&&(rows[i]!=5)&&(rows[i]!=10)&&(rows[i]!=11))
815  BOOST_CHECK_EQUAL(vals[i],-2.);
816 
817  sys->solution()->debug_data(vals);
818  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,-3.);
819 
820  sys->rhs()->debug_data(vals);
821  for (int i=0; i<vals.size(); i++)
822  {
823  if ((i==10)||(i==11)) { BOOST_CHECK_EQUAL(vals[i],0.); }
824  else { BOOST_CHECK_EQUAL(vals[i],-4.); }
825  }
826  }
827 
828  // diagonal access check
829  sys->reset();
830  std::vector<Real> diag(blockcol_size*neq);
831  for (int i=0; i<diag.size(); i++) diag[i]=i;
832  sys->set_diagonal(diag);
833  sys->add_diagonal(diag);
834  diag.clear();
835  sys->get_diagonal(diag);
836  BOOST_CHECK_EQUAL(diag.size(),blockcol_size*neq);
837  sys->matrix()->debug_data(rows,cols,vals);
838  for (int i=0; i<(const int)vals.size(); i++)
839  {
840  if (cp.isUpdatable()[rows[i]/neq]) { BOOST_CHECK_EQUAL(diag[rows[i]],2.*rows[i]); }
841  else { BOOST_CHECK_EQUAL(diag[rows[i]],0.); }
842  if (rows[i]==cols[i]) { BOOST_CHECK_EQUAL(vals[i],diag[rows[i]]); }
843  else { BOOST_CHECK_EQUAL(vals[i],0.); }
844  }
845 
846  // test swapping rhs and sol
847  boost::shared_ptr<LSS::System> sys2(common::allocate_component<LSS::System>("sys2"));
848  sys->options().option("matrix_builder").change_value(matrix_builder);
849  build_system(*sys2,cp);
850  BOOST_CHECK_EQUAL(sys2->is_created(),true);
851  BOOST_CHECK_EQUAL(sys2->solvertype(),solvertype);
852  sys->reset(1.);
853  sys2->reset(2.);
854  boost::shared_ptr<Matrix> sw_mat = boost::dynamic_pointer_cast<Matrix>(sys->remove_component("Matrix"));
855  boost::shared_ptr<Vector> sw_sol = boost::dynamic_pointer_cast<Vector>(sys2->remove_component("Solution"));
856  boost::shared_ptr<Vector> sw_rhs = boost::dynamic_pointer_cast<Vector>(sys2->remove_component("RHS"));
857  sys->swap(sw_mat,sw_sol,sw_rhs);
858  sys->matrix()->debug_data(rows,cols,vals);
859  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,1.);
860  sys->solution()->debug_data(vals);
861  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,2.);
862  sys->rhs()->debug_data(vals);
863  BOOST_FOREACH(Real i, vals) BOOST_CHECK_EQUAL(i,2.);
864 
865  // check destroy
866  sys->destroy();
867  BOOST_CHECK_EQUAL(sys->is_created(),false);
868 }
869 
871 
872 BOOST_AUTO_TEST_CASE( solve_system )
873 {
874 
875 // THE SERIAL IS EQUIVALENT WITH THE FOLLOWING OCTAVE/MATLAB CODE
876 // A = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
877 // b = [1; 1; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 10; 10]
878 // inv(A)*b
879 // WHICH RESULTS IN GID ORDER:
880 
881  std::vector<Real> refvals(0);
882  refvals +=
883  1.00000000000000e+00,
884  1.00000000000000e+00,
885  -1.35789473684210e+01,
886  -1.35789473684210e+01,
887  -7.78947368421052e+00,
888  -7.78947368421052e+00,
889  9.68421052631579e+00,
890  9.68421052631579e+00,
891  1.26315789473684e+01,
892  1.26315789473684e+01,
893  -3.36842105263158e+00,
894  -3.36842105263158e+00,
895  -1.43157894736842e+01,
896  -1.43157894736842e+01,
897  -3.78947368421053e+00,
898  -3.78947368421052e+00,
899  1.24210526315789e+01,
900  1.24210526315789e+01,
901  1.00000000000000e+01,
902  1.00000000000000e+01;
903 
904  // commpattern
905  if (irank==0)
906  {
907  gid += 0,1,2,3,4;
908  rank_updatable += 0,0,0,0,1;
909  } else {
910  gid += 3,4,5,6,7,8,9;
911  rank_updatable += 0,1,1,1,1,1,1;
912  }
913  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
914  common::PE::CommPattern& cp = *cp_ptr;
915  cp.insert("gid",gid,1,false);
916  cp.setup(Handle<common::PE::CommWrapper>(cp.get_child("gid")),rank_updatable);
917 
918  // lss
919  if (irank==0)
920  {
921  node_connectivity += 0,1,0,1,2,1,2,3,2,3,4,3,4;
922  starting_indices += 0,2,5,8,11,13;
923  } else {
924  node_connectivity += 0,1,0,1,2,1,2,3,2,3,4,3,4,5,4,5,6,5,6;
925  starting_indices += 0,2,5,8,11,14,17,19;
926  }
927  boost::shared_ptr<System> sys(common::allocate_component<System>("sys"));
928  sys->options().option("matrix_builder").change_value(matrix_builder);
929  sys->create(cp,2,node_connectivity,starting_indices);
930 
931  sys->solution_strategy()->options().set("compute_residual", true);
932  sys->solution_strategy()->options().set("verbosity_level", 3);
933  sys->solution_strategy()->access_component("Parameters")->options().set("preconditioner_type", std::string("None"));
934  sys->solution_strategy()->access_component("Parameters/LinearSolverTypes/Belos/SolverTypes/BlockGMRES")->options().set("verbosity", 1);
935 
936  // set intital values and boundary conditions
937  sys->matrix()->reset(-0.5);
938  sys->solution()->reset(1.);
939  sys->rhs()->reset(0.);
940  if (irank==0)
941  {
942  std::vector<Real> diag(10,1.);
943  sys->set_diagonal(diag);
944  sys->dirichlet(0,0,1.);
945  sys->dirichlet(0,1,1.);
946  } else {
947  std::vector<Real> diag(14,1.);
948  sys->set_diagonal(diag);
949  sys->dirichlet(6,0,10.);
950  sys->dirichlet(6,1,10.);
951  }
952 
953  // solve and check
954  sys->solve();
955  std::vector<Real> vals;
956  sys->solution()->debug_data(vals);
957  for (int i=0; i<vals.size(); i++)
958  if (cp.isUpdatable()[i/neq])
959  BOOST_CHECK_CLOSE( vals[i], refvals[gid[i/neq]*neq], 1e-8);
960 
961 }
962 
964 
965 BOOST_AUTO_TEST_CASE( solve_system_blocked )
966 {
967 
968 // THE SERIAL IS EQUIVALENT WITH THE FOLLOWING OCTAVE/MATLAB CODE
969 // A = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, 1, -0.5, -0.5, -0.5; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5, -0.5, -0.5, 1, -0.5, -0.5; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
970 // b = [1; 1; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 10; 10]
971 // inv(A)*b
972 // WHICH RESULTS IN GID ORDER:
973 
974  try
975  {
976  std::vector<Real> refvals(0);
977  refvals +=
978  1.00000000000000e+00,
979  1.00000000000000e+00,
980  -1.35789473684210e+01,
981  -1.35789473684210e+01,
982  -7.78947368421052e+00,
983  -7.78947368421052e+00,
984  9.68421052631579e+00,
985  9.68421052631579e+00,
986  1.26315789473684e+01,
987  1.26315789473684e+01,
988  -3.36842105263158e+00,
989  -3.36842105263158e+00,
990  -1.43157894736842e+01,
991  -1.43157894736842e+01,
992  -3.78947368421053e+00,
993  -3.78947368421052e+00,
994  1.24210526315789e+01,
995  1.24210526315789e+01,
996  1.00000000000000e+01,
997  1.00000000000000e+01;
998 
999  // commpattern
1000  if (irank==0)
1001  {
1002  gid += 0,1,2,3,4;
1003  rank_updatable += 0,0,0,0,1;
1004  } else {
1005  gid += 3,4,5,6,7,8,9;
1006  rank_updatable += 0,1,1,1,1,1,1;
1007  }
1008  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
1009  common::PE::CommPattern& cp = *cp_ptr;
1010  cp.insert("gid",gid,1,false);
1011  cp.setup(Handle<common::PE::CommWrapper>(cp.get_child("gid")),rank_updatable);
1012 
1013  // lss
1014  if (irank==0)
1015  {
1016  node_connectivity += 0,1,0,1,2,1,2,3,2,3,4,3,4;
1017  starting_indices += 0,2,5,8,11,13;
1018  } else {
1019  node_connectivity += 0,1,0,1,2,1,2,3,2,3,4,3,4,5,4,5,6,5,6;
1020  starting_indices += 0,2,5,8,11,14,17,19;
1021  }
1022  boost::shared_ptr<System> sys(common::allocate_component<System>("sys"));
1023  sys->options().option("matrix_builder").change_value(matrix_builder);
1024  boost::shared_ptr<math::VariablesDescriptor> vars = common::allocate_component<math::VariablesDescriptor>("vars");
1025  vars->options().set("dimension", 1u);
1026 
1027  vars->push_back("var1", cf3::math::VariablesDescriptor::Dimensionalities::SCALAR);
1028  vars->push_back("var2", cf3::math::VariablesDescriptor::Dimensionalities::SCALAR);
1029  sys->create_blocked(cp,*vars,node_connectivity,starting_indices);
1030 
1031  sys->solution_strategy()->options().set("compute_residual", true);
1032  sys->solution_strategy()->options().set("verbosity_level", 3);
1033  sys->solution_strategy()->access_component("Parameters")->options().set("preconditioner_type", std::string("None"));
1034  sys->solution_strategy()->access_component("Parameters/LinearSolverTypes/Belos/SolverTypes/BlockGMRES")->options().set("verbosity", 1);
1035  sys->solution_strategy()->access_component("Parameters/LinearSolverTypes/Belos/SolverTypes/BlockGMRES")->options().set("convergence_tolerance", 0.1);
1036 
1037  // set intital values and boundary conditions
1038  sys->matrix()->reset(-0.5);
1039  sys->solution()->reset(1.);
1040  sys->rhs()->reset(0.);
1041  if (irank==0)
1042  {
1043  std::vector<Real> diag(10,1.);
1044  sys->set_diagonal(diag);
1045  sys->dirichlet(0,0,1.);
1046  sys->dirichlet(0,1,1.);
1047  } else {
1048  std::vector<Real> diag(14,1.);
1049  sys->set_diagonal(diag);
1050  sys->dirichlet(6,0,10.);
1051  sys->dirichlet(6,1,10.);
1052  }
1053 
1054  // solve and check
1055  sys->solve();
1056  std::vector<Real> vals;
1057  sys->solution()->debug_data(vals);
1058  for (int i=0; i<vals.size(); i++)
1059  if (cp.isUpdatable()[i/neq])
1060  BOOST_CHECK_CLOSE( vals[i], refvals[gid[i/neq]*neq], 1e-8);
1061  }
1062  catch(common::NotImplemented&)
1063  {
1064  CFinfo << "skipping blocked solve" << CFendl;
1065  }
1066 }
1067 
1069 
1070 BOOST_AUTO_TEST_CASE( solve_system_laplacian )
1071 {
1072  // commpattern
1073  if (irank==0)
1074  {
1075  gid += 0,1,2,3;
1076  rank_updatable += 0,0,0,1;
1077  } else {
1078  gid += 2,3,4,5,6;
1079  rank_updatable += 0,1,1,1,1;
1080  }
1081  boost::shared_ptr<common::PE::CommPattern> cp_ptr = common::allocate_component<common::PE::CommPattern>("commpattern");
1082  common::PE::CommPattern& cp = *cp_ptr;
1083  cp.insert("gid",gid,1,false);
1084  cp.setup(Handle<common::PE::CommWrapper>(cp.get_child("gid")),rank_updatable);
1085 
1086  // lss
1087  if (irank==0)
1088  {
1089  node_connectivity += 0,1,0,1,2,1,2,3,2,3;
1090  starting_indices += 0,2,5,8,10;
1091  } else {
1092  node_connectivity += 0,1,0,1,2,1,2,3,2,3,4,3,4;
1093  starting_indices += 0,2,5,8,11,13;
1094  }
1095  boost::shared_ptr<System> sys(common::allocate_component<System>("sys"));
1096  sys->options().option("matrix_builder").change_value(matrix_builder);
1097  sys->create(cp,1,node_connectivity,starting_indices);
1098 
1099  sys->solution_strategy()->options().set("compute_residual", true);
1100  sys->solution_strategy()->options().set("verbosity_level", 3);
1101  sys->solution_strategy()->options().set("print_settings", false);
1102  sys->solution_strategy()->access_component("Parameters")->options().set("preconditioner_type", std::string("None"));
1103  sys->solution_strategy()->access_component("Parameters")->options().set("linear_solver_type", std::string("Amesos"));
1104  //sys->solution_strategy()->access_component("Parameters/LinearSolverTypes/Belos/SolverTypes/BlockGMRES")->options().set("verbosity", 1);
1105 
1106  // set intital values and boundary conditions
1107  sys->matrix()->reset(1.);
1108  sys->solution()->reset(0.);
1109  sys->rhs()->reset(0.);
1110  if (irank==0)
1111  {
1112  std::vector<Real> diag(4,-2.);
1113  sys->set_diagonal(diag);
1114  sys->dirichlet(0,0,10.);
1115  } else {
1116  std::vector<Real> diag(5,-2.);
1117  sys->set_diagonal(diag);
1118  sys->dirichlet(4,0,16.);
1119  }
1120 
1121  // solve and check
1122  sys->solve();
1123  sys->solution()->print_native(std::cout);
1124 // std::vector<Real> vals;
1125 // sys->solution()->debug_data(vals);
1126 // for (int i=0; i<vals.size(); i++)
1127 // if (cp.isUpdatable()[i/neq])
1128 // BOOST_CHECK_CLOSE( vals[i], refvals[gid[i/neq]*neq], 1e-8);
1129 
1130 }
1131 
1133 
1134 BOOST_AUTO_TEST_CASE( finalize_mpi )
1135 {
1136  CFinfo.setFilterRankZero(true);
1137  common::PE::Comm::instance().finalize();
1138  BOOST_CHECK_EQUAL(common::PE::Comm::instance().is_active(),false);
1139 }
1140 
1142 
1143 BOOST_AUTO_TEST_SUITE_END()
1144 
1145 
cf3::math::LSS::Vector::debug_data
virtual void debug_data(std::vector< Real > &values)=0
LSSAtomicFixture::build_system
void build_system(LSS::System &sys, common::PE::CommPattern &cp)
build a test system
Definition: utest-lss-atomic.cpp:87
CFinfo
#define CFinfo
these are always defined
Definition: Log.hpp:104
cf3::math::LSS::System::create
void create(cf3::common::PE::CommPattern &cp, Uint neq, std::vector< Uint > &node_connectivity, std::vector< Uint > &starting_indices, const std::vector< Uint > &periodic_links_nodes=std::vector< Uint >(), const std::vector< bool > &periodic_links_active=std::vector< bool >())
Definition: System.cpp:82
cf3::math::LSS::Vector::add_rhs_values
virtual void add_rhs_values(const BlockAccumulator &values)=0
Add a list of values to rhs.
cf3::math::LSS::BlockAccumulator::resize
void resize(Uint numnodes, Uint numeqs)
setting up sizes
Definition: BlockAccumulator.hpp:32
cf3::math::LSS::Vector::set_rhs_values
virtual void set_rhs_values(const BlockAccumulator &values)=0
Set a list of values to rhs.
cf3::math::LSS::Vector::destroy
virtual void destroy()=0
Deallocate underlying data.
cf3::math::LSS::Matrix::reset
virtual void reset(Real reset_to=0.)=0
Reset Matrix.
cf3::Handle
Safe pointer to an object. This is the supported method for referring to components.
Definition: Handle.hpp:39
System.hpp
This header collects all the headers needed for the linear system solver, also including configure-ti...
cf3::math
Basic Classes for Mathematical applications used by COOLFluiD.
Definition: AnalyticalFunction.cpp:24
Environment.hpp
cf3::math::LSS::Matrix::get_values
virtual void get_values(BlockAccumulator &values)=0
Add a list of values.
cf3::common::PE::CommPattern::setup
void setup(const Handle< CommWrapper > &gid, std::vector< Uint > &rank)
cf3::math::LSS::Matrix::set_values
virtual void set_values(const BlockAccumulator &values)=0
Set a list of values.
cf3::math::LSS::Matrix::blockrow_size
virtual const Uint blockrow_size()=0
Accessor to the number of block rows.
cf3::math::LSS::BlockAccumulator::mat
Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > mat
Definition: BlockAccumulator.hpp:61
LSSAtomicFixture::build_commpattern
void build_commpattern(common::PE::CommPattern &cp)
create a test commpattern
Definition: utest-lss-atomic.cpp:72
cf3::math::LSS::Vector::blockrow_size
virtual const Uint blockrow_size()=0
Accessor to the number of block rows.
cf3::math::LSS::Matrix::destroy
virtual void destroy()=0
Deallocate underlying data.
cf3::math::LSS::Matrix::set_value
virtual void set_value(const Uint icol, const Uint irow, const Real value)=0
Set value at given location in the matrix.
cf3::math::LSS::Matrix::get_diagonal
virtual void get_diagonal(std::vector< Real > &diag)=0
Get the diagonal.
CFendl
#define CFendl
Definition: Log.hpp:109
cf3::common::PE::CommPattern::insert
void insert(const std::string &name, T *&data, const int size, const unsigned int stride=1, const bool needs_update=true)
Definition: CommPattern.hpp:139
cf3::math::LSS::Matrix::is_created
virtual const bool is_created()=0
Accessor to the state of create.
cf3::physics::Consts::e
Real e()
Definition of the Unit charge [C].
Definition: Consts.hpp:30
cf3::math::LSS::Matrix::add_values
virtual void add_values(const BlockAccumulator &values)=0
cf3::math::LSS::Vector::get_sol_values
virtual void get_sol_values(BlockAccumulator &values)=0
Get a list of values from sol.
LSSAtomicFixture::rank_updatable
std::vector< Uint > rank_updatable
Definition: utest-lss-atomic.cpp:119
cf3::math::LSS::Matrix::set_diagonal
virtual void set_diagonal(const std::vector< Real > &diag)=0
Set the diagonal.
cf3::math::LSS::BlockAccumulator::indices
std::vector< Uint > indices
local numbering of the unknowns
Definition: BlockAccumulator.hpp:74
cf3::math::LSS::Vector::set_sol_values
virtual void set_sol_values(const BlockAccumulator &values)=0
Set a list of values to sol.
cf3::math::LSS::Matrix::set_row
virtual void set_row(const Uint iblockrow, const Uint ieq, Real diagval, Real offdiagval)=0
Set a row, diagonal and off-diagonals values separately (dirichlet-type boundaries) ...
cf3::math::LSS::Vector::get_rhs_values
virtual void get_rhs_values(BlockAccumulator &values)=0
Get a list of values from rhs.
VariablesDescriptor.hpp
cf3::math::LSS::Matrix::tie_blockrow_pairs
virtual void tie_blockrow_pairs(const Uint iblockrow_to, const Uint iblockrow_from)=0
Add one line to another and tie to it via dirichlet-style (applying periodicity)
cf3::math::LSS::Matrix::get_column_and_replace_to_zero
virtual void get_column_and_replace_to_zero(const Uint iblockcol, Uint ieq, std::vector< Real > &values)=0
cf3::math::LSS::Vector::set_value
virtual void set_value(const Uint irow, const Real value)=0
Set value at given location in the matrix.
cf3::common::Component::get_child
Handle< Component > get_child(const std::string &name)
Definition: Component.cpp:441
cf3
Top-level namespace for coolfluid.
Definition: Action.cpp:18
LSSAtomicFixture::~LSSAtomicFixture
~LSSAtomicFixture()
common tear-down for each test case
Definition: utest-lss-atomic.cpp:67
cf3::math::LSS::Matrix::add_diagonal
virtual void add_diagonal(const std::vector< Real > &diag)=0
Add to the diagonal.
cf3::math::LSS::Matrix::add_value
virtual void add_value(const Uint icol, const Uint irow, const Real value)=0
Add value at given location in the matrix.
cf3::math::LSS::Matrix::blockcol_size
virtual const Uint blockcol_size()=0
Accessor to the number of block columns.
cf3::math::LSS::Vector::get_value
virtual void get_value(const Uint irow, Real &value)=0
Get value at given location in the matrix.
cf3::math::LSS::BlockAccumulator::sol
RealVector sol
accessor to blockaccumulator's solution vector
Definition: BlockAccumulator.hpp:68
LSSAtomicFixture::starting_indices
std::vector< Uint > starting_indices
Definition: utest-lss-atomic.cpp:123
cf3::math::LSS::Vector::is_created
virtual const bool is_created()=0
Accessor to the state of create.
LSSAtomicFixture::solvertype
std::string solvertype
main solver selector
Definition: utest-lss-atomic.cpp:107
LSSAtomicFixture::irank
int irank
constructor builds
Definition: utest-lss-atomic.cpp:111
cf3::math::LSS::Matrix::debug_data
virtual void debug_data(std::vector< Uint > &row_indices, std::vector< Uint > &col_indices, std::vector< Real > &values)=0
cf3::math::LSS::BlockAccumulator::reset
void reset(Real reset_to=0.)
reset the values to the value of reset_to
Definition: BlockAccumulator.hpp:42
Log.hpp
cf3::math::LSS::Matrix::neq
virtual const Uint neq()=0
Accessor to the number of equations.
cf3::math::LSS::Vector::print
virtual void print(common::LogStream &stream)=0
Print to wherever.
LSSAtomicFixture::LSSAtomicFixture
LSSAtomicFixture()
common setup for each test case
Definition: utest-lss-atomic.cpp:40
LSSAtomicFixture::gid
std::vector< Uint > gid
commpattern builds
Definition: utest-lss-atomic.cpp:118
LSSAtomicFixture::matrix_builder
std::string matrix_builder
Definition: utest-lss-atomic.cpp:108
BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(init_mpi)
Definition: utest-lss-atomic.cpp:135
LSSAtomicFixture::node_connectivity
std::vector< Uint > node_connectivity
system builds
Definition: utest-lss-atomic.cpp:122
cf3::math::LSS::Matrix::print
virtual void print(common::LogStream &stream)=0
Print to wherever.
cf3::math::LSS::Vector::add_value
virtual void add_value(const Uint irow, const Real value)=0
Add value at given location in the matrix.
cf3::math::LSS::Vector::solvertype
virtual const std::string solvertype()=0
Accessor to solver type.
cf3::math::LSS::Vector::reset
virtual void reset(Real reset_to=0.)=0
Reset Vector.
cf3::math::LSS::BlockAccumulator::rhs
RealVector rhs
accessor to blockaccumulator's right hand side vector
Definition: BlockAccumulator.hpp:71
cf3::math::LSS::Matrix::solvertype
virtual const std::string solvertype()=0
Accessor to solver type.
cf3::math::LSS::Vector::add_sol_values
virtual void add_sol_values(const BlockAccumulator &values)=0
Add a list of values to sol.
cf3::math::LSS::Vector::neq
virtual const Uint neq()=0
Accessor to the number of equations.
cf3::math::LSS::Matrix::get_value
virtual void get_value(const Uint icol, const Uint irow, Real &value)=0
Get value at given location in the matrix.
FromHere
#define FromHere()
Definition: CodeLocation.hpp:45
Generated on Sun Jun 14 2015 21:20:16 for COOLFluiD by doxygen 1.8.9.1
Send comments to:
COOLFluiD Web Admin