MOOSE - Multiscale Object Oriented Simulation Environment
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
testDiffusion.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 ** This program is part of 'MOOSE', the
3 ** Messaging Object Oriented Simulation Environment.
4 ** Copyright (C) 2003-2014 Upinder S. Bhalla. and NCBS
5 ** It is made available under the terms of the
6 ** GNU Lesser General Public License version 2.1
7 ** See the file COPYING.LIB for the full notice.
8 **********************************************************************/
9 /*
10 #include <vector>
11 #include <algorithm>
12 #include <cassert>
13 #include <functional>
14 #include <iostream>
15 #include <iomanip>
16 //#include "/usr/include/gsl/gsl_linalg.h"
17 using namespace std;
18 */
19 #ifdef USE_GSL
20 #include <gsl/gsl_linalg.h>
21 #endif
22 #include "header.h"
23 #include "../basecode/SparseMatrix.h"
24 #include "FastMatrixElim.h"
25 #include "../shell/Shell.h"
26 
27 
28 
29 double checkAns(
30  const double* m, unsigned int numCompts,
31  const double* ans, const double* rhs )
32 {
33  vector< double > check( numCompts, 0.0 );
34  for ( unsigned int i = 0; i < numCompts; ++i ) {
35  for ( unsigned int j = 0; j < numCompts; ++j )
36  check[i] += m[i*numCompts + j] * ans[j];
37  }
38  double ret = 0.0;
39  for ( unsigned int i = 0; i < numCompts; ++i )
40  ret += (check[i] - rhs[i]) * (check[i] - rhs[i] );
41  return ret;
42 }
43 
44 
45 
46 void testFastMatrixElim()
47 {
48 
49 
50 /*
51 2 11
52  1 4
53  3 10
54  9 5
55  6
56  7
57  8
58 
59  1 2 3 4 5 6 7 8 9 10 11
60 1 x x x x
61 2 x x
62 3 x x x x
63 4 x x x x
64 5 x x x x
65 6 x x x x
66 7 x x x
67 8 x x
68 9 x x x x
69 10 x x
70 11 x x
71  static double test[] = {
72  1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0,
73  5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0,
74  7, 0, 8, 9, 0, 0, 0, 0, 10, 0, 0,
75  11, 0, 12, 13, 0, 0, 0, 0, 0, 0, 14,
76  0, 0, 0, 0, 15, 16, 0, 0, 17, 18, 0,
77  0, 0, 0, 0, 19, 20, 21, 0, 22, 0, 0,
78  0, 0, 0, 0, 0, 23, 24, 25, 0, 0, 0,
79  0, 0, 0, 0, 0, 0, 26, 27, 0, 0, 0,
80  0, 0, 28, 0, 29, 30, 0, 0, 31, 0, 0,
81  0, 0, 0, 0, 32, 0, 0, 0, 0, 33, 0,
82  0, 0, 0, 34, 0, 0, 0, 0, 0, 0, 35,
83  };
84  const unsigned int numCompts = 11;
85 // static unsigned int parents[] = { 3,1,9,3,6,7,8,-1,6,5,4 };
86  static unsigned int parents[] = { 2,0,8,2,5,6,7,-1,5,4,3 };
87 */
88 
89 /*
90 1 3
91  2 4
92  7 5
93  8 6
94  9
95  10
96  11
97 
98  1 2 3 4 5 6 7 8 9 10 11
99 1 x x
100 2 x x x
101 3 x x
102 4 x x x
103 5 x x
104 6 x x x
105 7 x x x x
106 8 x x x
107 9 x x x x
108 10 x x x
109 11 x x
110  static double test[] = {
111  1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
112  3, 4, 0, 0, 0, 0, 5, 0, 0, 0, 0,
113  0, 0, 6, 7, 0, 0, 0, 0, 0, 0, 0,
114  0, 0, 8, 9, 0, 0, 10, 0, 0, 0, 0,
115  0, 0, 0, 0, 11, 12, 0, 0, 0, 0, 0,
116  0, 0, 0, 0, 13, 14, 0, 0, 15, 0, 0,
117  0, 16, 0, 17, 0, 0, 18, 19, 0, 0, 0,
118  0, 0, 0, 0, 0, 0, 20, 21, 22, 0, 0,
119  0, 0, 0, 0, 0, 23, 0, 24, 25, 26, 0,
120  0, 0, 0, 0, 0, 0, 0, 0, 27, 28, 29,
121  0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 31,
122  };
123  const unsigned int numCompts = 11;
124  static unsigned int parents[] = { 1,6,3,6,5,8,7,8,9,10,-1};
125 */
126 
127 /*
128 Linear cable, 12 segments.
129 */
130 
131  static double test[] = {
132  1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
133  3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
134  0, 6, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0,
135  0, 0, 9, 10, 11, 0, 0, 0, 0, 0, 0, 0,
136  0, 0, 0, 12, 13, 14, 0, 0, 0, 0, 0, 0,
137  0, 0, 0, 0, 15, 16, 17, 0, 0, 0, 0, 0,
138  0, 0, 0, 0, 0, 18, 19, 20, 0, 0, 0, 0,
139  0, 0, 0, 0, 0, 0, 21, 22, 23, 0, 0, 0,
140  0, 0, 0, 0, 0, 0, 0, 24, 25, 26, 0, 0,
141  0, 0, 0, 0, 0, 0, 0, 0, 27, 28, 29, 0,
142  0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 31, 32,
143  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 33, 34,
144  };
145  const unsigned int numCompts = 12;
146  static unsigned int parents[] = { 1,2,3,4,5,6,7,8,9,10,11, unsigned(-1)};
147 
148  /*
149  static double test[] = {
150  1, 2,
151  3, 4
152  };
153  const unsigned int numCompts = 2;
154  static double test[] = {
155  1, 2, 0, 0,
156  3, 4, 5, 0,
157  0, 6, 7, 8,
158  0, 0, 9, 10
159  };
160  const unsigned int numCompts = 4;
161  static double test[] = {
162  1, 2, 0, 0, 0, 0,
163  3, 4, 5, 0, 0, 0,
164  0, 6, 7, 8, 0, 0,
165  0, 0, 9, 10, 11, 0,
166  0, 0, 0, 12, 13, 14,
167  0, 0, 0, 0, 15, 16,
168  };
169  const unsigned int numCompts = 6;
170  static double test[] = {
171  1, 2, 0, 0, 0, 0,
172  3, 4, 0, 0, 1, 0,
173  0, 0, 7, 8, 0, 0,
174  0, 0, 9, 10, 11, 0,
175  0, 1, 0, 12, 13, 14,
176  0, 0, 0, 0, 15, 16,
177  };
178  const unsigned int numCompts = 6;
179  */
180  // testSorting(); // seems to work fine.
181  FastMatrixElim fe;
182  vector< Triplet< double > > fops;
183  fe.makeTestMatrix( test, numCompts );
184  // fe.print();
185  vector< unsigned int > parentVoxel;
186  vector< unsigned int > lookupOldRowsFromNew;
187  parentVoxel.insert( parentVoxel.begin(), &parents[0], &parents[numCompts] );
188  fe.hinesReorder( parentVoxel, lookupOldRowsFromNew );
189  /*
190  */
191  /*
192  vector< unsigned int > shuf;
193  for ( unsigned int i = 0; i < numCompts; ++i )
194  shuf.push_back( i );
195  shuf[0] = 1;
196  shuf[1] = 0;
197  fe.shuffleRows( shuf );
198  */
199  // fe.print();
200  FastMatrixElim foo = fe;
201 
202  vector< unsigned int > diag;
203  vector< double > diagVal;
204  fe.buildForwardElim( diag, fops );
205  // fe.print();
206  fe.buildBackwardSub( diag, fops, diagVal );
207  vector< double > y( numCompts, 1.0 );
208  vector< double > ones( numCompts, 1.0 );
209  FastMatrixElim::advance( y, fops, diagVal );
210  /*
211  for ( unsigned int i = 0; i < numCompts; ++i )
212  cout << "y" << i << "]= " << y[i] << endl;
213  */
214 
215  // Here we verify the answer
216 
217  vector< double > alle;
218  for( unsigned int i = 0; i < numCompts; ++i ) {
219  for( unsigned int j = 0; j < numCompts; ++j ) {
220  alle.push_back( foo.get( i, j ) );
221  }
222  }
223  // cout << "myCode: " << checkAns( &alle[0], numCompts, &y[0], &ones[0] ) << endl;
224 
225  assert( checkAns( &alle[0], numCompts, &y[0], &ones[0] ) < 1e-25 );
226 
227 #if USE_GSL
228  // Here we do the gsl test.
230  vector< double > temp( &test[0], &test[numCompts*numCompts] );
231  gsl_matrix_view m = gsl_matrix_view_array( &temp[0], numCompts, numCompts );
232 
233  vector< double > z( numCompts, 1.0 );
234  gsl_vector_view b = gsl_vector_view_array( &z[0], numCompts );
235  gsl_vector* x = gsl_vector_alloc( numCompts );
236  int s;
237  gsl_permutation* p = gsl_permutation_alloc( numCompts );
238  gsl_linalg_LU_decomp( &m.matrix, p, &s );
239  gsl_linalg_LU_solve( &m.matrix, p, &b.vector, x);
240  vector< double > gslAns( numCompts );
241  for ( unsigned int i = 0; i < numCompts; ++i ) {
242  gslAns[i] = gsl_vector_get( x, i );
243  // cout << "x[" << i << "]= " << gslAns[i] << endl;
244  }
245  // cout << "GSL: " << checkAns( test, numCompts, &gslAns[0], &ones[0] ) << endl;
246  assert( checkAns( test, numCompts, &gslAns[0], &ones[0] ) < 1e-25 );
247  gsl_permutation_free( p );
248  gsl_vector_free( x );
249  cout << "." << flush;
250 #endif
251 }
252 
253 void testSorting()
254 {
255  static unsigned int k[] = {20,40,60,80,100,10,30,50,70,90};
256  static double d[] = {1,2,3,4,5,6,7,8,9,10};
257  vector< unsigned int > col;
258  col.insert( col.begin(), k, k+10);
259  vector< double > entry;
260  entry.insert( entry.begin(), d, d+10);
261  sortByColumn( col, entry );
262 
263  for ( unsigned int i = 0; i < col.size(); ++i )
264  assert( col[i] == (i + 1) * 10 );
265 
266  assert( entry[0] == 6 );
267  assert( entry[1] == 1 );
268  assert( entry[2] == 7 );
269  assert( entry[3] == 2 );
270  assert( entry[4] == 8 );
271  assert( entry[5] == 3 );
272  assert( entry[6] == 9 );
273  assert( entry[7] == 4 );
274  assert( entry[8] == 10 );
275  assert( entry[9] == 5 );
276 
277  /*
278  cout << "testing sorting\n";
279  for ( unsigned int i = 0; i < col.size(); ++i ) {
280  cout << "d[" << i << "]= " << k[i] <<
281  ", col[" << i << "]= " << col[i] << ", e=" << entry[i] << endl;
282  }
283  cout << endl;
284  */
285  cout << "." << flush;
286 }
287 
288 void testSetDiffusionAndTransport()
289 {
290  static double test[] = {
291  0, 2, 0, 0, 0, 0,
292  1, 0, 2, 0, 0, 0,
293  0, 1, 0, 2, 0, 0,
294  0, 0, 1, 0, 2, 0,
295  0, 0, 0, 1, 0, 2,
296  0, 0, 0, 0, 1, 0,
297  };
298  const unsigned int numCompts = 6;
299  FastMatrixElim fm;
300  fm.makeTestMatrix( test, numCompts );
301  vector< unsigned int > parentVoxel( numCompts );
302  parentVoxel[0] = -1;
303  parentVoxel[1] = 0;
304  parentVoxel[2] = 1;
305  parentVoxel[3] = 2;
306  parentVoxel[4] = 3;
307  parentVoxel[5] = 4;
308 
309  // cout << endl;
310  // fm.print();
311  // cout << endl;
312  // fm.printInternal();
313  fm.setDiffusionAndTransport( parentVoxel, 1, 10, 0.1 );
314  // cout << endl;
315  // fm.print();
316  // cout << endl;
317  // fm.printInternal();
318 
319  for( unsigned int i =0; i < numCompts; ++i ) {
320  unsigned int start = 0;
321  if ( i > 0 )
322  start = i - 1;
323  for( unsigned int j = start ; j < i+1 && j < numCompts ; ++j ) {
324  if ( i == j + 1 )
325  assert( doubleEq( fm.get( i, j ), 0.1 ) );
326  else if ( i + 1 == j ) {
327  assert( doubleEq( fm.get( i, j ), 2.2 ) );
328  } else if ( i == j ) {
329  if ( i == 0 )
330  assert( doubleEq( fm.get( i, j ), 0.8 ) );
331  else if ( i == numCompts - 1 )
332  assert( doubleEq( fm.get( i, j ), -0.1 ) );
333  else
334  assert( doubleEq( fm.get( i, j ), -0.3 ) );
335  }
336  }
337  }
338  cout << "." << flush;
339 }
340 
341 void testCylDiffn()
342 {
343  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
344  double len = 25e-6;
345  double r0 = 1e-6;
346  double r1 = 1e-6;
347  double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
348  double runtime = 10.0;
349  double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
350  double diffConst = 1.0e-12;
351  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
352  Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
353  Field< double >::set( cyl, "r0", r0 );
354  Field< double >::set( cyl, "r1", r1 );
355  Field< double >::set( cyl, "x0", 0 );
356  Field< double >::set( cyl, "x1", len );
357  Field< double >::set( cyl, "diffLength", diffLength );
358  unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
359  assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
360  Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
361  Field< double >::set( pool, "diffConst", diffConst );
362 
363  Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
364  Field< Id >::set( dsolve, "compartment", cyl );
365  // Next: build by doing reinit
366  s->doUseClock( "/model/dsolve", "process", 1 );
367  s->doSetClock( 1, dt );
368  Field< string >::set( dsolve, "path", "/model/cyl/pool" );
369  // Then find a way to test it.
370  vector< double > poolVec;
371  Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );
372  Field< double >::getVec( pool, "nInit", poolVec );
373  assert( poolVec.size() == ndc );
374  assert( doubleEq( poolVec[0], 1.0 ) );
375  assert( doubleEq( poolVec[1], 0.0 ) );
376 
377  vector< double > nvec =
378  LookupField< unsigned int, vector< double > >::get(
379  dsolve, "nVec", 0);
380  assert( nvec.size() == ndc );
381 
382  s->doReinit();
383  s->doStart( runtime );
384 
385  nvec = LookupField< unsigned int, vector< double > >::get(
386  dsolve, "nVec", 0);
387  Field< double >::getVec( pool, "n", poolVec );
388  assert( nvec.size() == poolVec.size() );
389  for ( unsigned int i = 0; i < nvec.size(); ++i )
390  assert( doubleEq( nvec[i], poolVec[i] ) );
391  /*
392  cout << endl;
393  for ( unsigned int i = 0; i < nvec.size(); ++i )
394  cout << nvec[i] << " ";
395  cout << endl;
396  */
397 
398  double dx = diffLength;
399  double err = 0.0;
400  double analyticTot = 0.0;
401  double myTot = 0.0;
402  for ( unsigned int i = 0; i < nvec.size(); ++i ) {
403  double x = i * dx + dx * 0.5;
404  // This part is the solution as a func of x,t.
405  double y = dx * // This part represents the init n of 1 in dx
406  ( 1.0 / sqrt( PI * diffConst * runtime ) ) *
407  exp( -x * x / ( 4 * diffConst * runtime ) );
408  err += ( y - nvec[i] ) * ( y - nvec[i] );
409  //cout << i << " " << x << " " << y << " " << conc[i] << endl;
410  analyticTot += y;
411  myTot += nvec[i];
412  }
413  assert( doubleEq( myTot, 1.0 ) );
414  // cout << "analyticTot= " << analyticTot << ", myTot= " << myTot << endl;
415  assert( err < 1.0e-5 );
416 
417 
418  s->doDelete( model );
419  cout << "." << flush;
420 }
421 
422 void testTaperingCylDiffn()
423 {
424  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
425  double len = 25e-6;
426  double r0 = 2e-6;
427  double r1 = 1e-6;
428  double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
429  double runtime = 10.0;
430  double dt = 0.1; // 0.2 is the highest dt for which the error is in bounds
431  double diffConst = 1.0e-12;
432  // Should set explicitly, currently during creation of DiffPoolVec
433  //double diffConst = 1.0e-12;
434  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
435  Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
436  Field< double >::set( cyl, "r0", r0 );
437  Field< double >::set( cyl, "r1", r1 );
438  Field< double >::set( cyl, "x0", 0 );
439  Field< double >::set( cyl, "x1", len );
440  Field< double >::set( cyl, "diffLength", diffLength );
441  unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
442  assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
443  Id pool = s->doCreate( "Pool", cyl, "pool", 1 );
444  Field< double >::set( pool, "diffConst", diffConst );
445 
446  Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
447  Field< Id >::set( dsolve, "compartment", cyl );
448  s->doUseClock( "/model/dsolve", "process", 1 );
449  s->doSetClock( 1, dt );
450  // Next: build by setting the path of the dsolve.
451  Field< string >::set( dsolve, "path", "/model/cyl/pool" );
452  // Then find a way to test it.
453  assert( pool.element()->numData() == ndc );
454  Field< double >::set( ObjId( pool, 0 ), "nInit", 1.0 );
455 
456  s->doReinit();
457  s->doStart( runtime );
458 
459  double myTot = 0.0;
460  vector< double > poolVec;
461  Field< double >::getVec( pool, "n", poolVec );
462  for ( unsigned int i = 0; i < poolVec.size(); ++i ) {
463  myTot += poolVec[i];
464  }
465  assert( doubleEq( myTot, 1.0 ) );
466 
467  s->doDelete( model );
468  cout << "." << flush;
469 }
470 
498 void testSmallCellDiffn()
499 {
500  Id makeCompt( Id parentCompt, Id parentObj,
501  string name, double len, double dia, double theta );
502  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
503  double len = 20e-6;
504  double dia = 10e-6;
505  double diffLength = 10e-6;
506  double dt = 1.0e-1;
507  double runtime = 100.0;
508  double diffConst = 1.0e-12;
509  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
510  Id soma = makeCompt( Id(), model, "soma", dia, dia, 90 );
511  Id dend = makeCompt( soma, model, "dend", len, 3e-6, 0 );
512  Id branch1 = makeCompt( dend, model, "branch1", len, 2e-6, 45.0 );
513  Id branch2 = makeCompt( dend, model, "branch2", len, 2e-6, -45.0 );
514  Id twig1 = makeCompt( branch1, model, "twig1", len, 1.5e-6, 90.0 );
515  Id twig2 = makeCompt( branch1, model, "twig2", len, 1.5e-6, 0.0 );
516 
517  Id nm = s->doCreate( "NeuroMesh", model, "neuromesh", 1 );
518  Field< double >::set( nm, "diffLength", diffLength );
519  Field< string >::set( nm, "geometryPolicy", "cylinder" );
520  Field< string >::set( nm, "subTreePath", "/model/#" );
521  //SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );
522  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
523  assert( ns == 6 );
524  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
525  assert( ndc == 11 );
526  Id pool1 = s->doCreate( "Pool", nm, "pool1", 1 );
527  Field< double >::set( pool1, "diffConst", diffConst );
528  Id pool2 = s->doCreate( "Pool", nm, "pool2", 1 );
529  Field< double >::set( pool2, "diffConst", diffConst );
530  Id pool3 = s->doCreate( "Pool", nm, "pool3", 1 );
531  Field< double >::set( pool3, "diffConst", 0.0 );
532 
533  Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
534  Field< Id >::set( dsolve, "compartment", nm );
535  s->doUseClock( "/model/dsolve", "process", 1 );
536  s->doSetClock( 1, dt );
537  // Next: build diffusion by setting path
538  Field< string >::set( dsolve, "path", "/model/neuromesh/pool#" );
539 
540  vector< double > nvec =
541  LookupField< unsigned int, vector< double > >::get(
542  dsolve, "nVec", 0);
543  assert( nvec.size() == ndc );
544  assert( pool1.element()->numData() == ndc );
545  Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
546  Field< double >::set( ObjId( pool2, ndc - 1 ), "nInit", 2.0 );
547  Field< double >::set( ObjId( pool3, 0 ), "nInit", 3.0 );
548 
549  s->doReinit();
550  nvec = LookupField< unsigned int, vector< double > >::get(
551  dsolve, "nVec", 0);
552  assert( doubleEq( nvec[0], 1.0 ) );
553  assert( doubleEq( nvec[1], 0.0 ) );
554  s->doStart( runtime );
555 
556  nvec = LookupField< unsigned int, vector< double > >::get(
557  dsolve, "nVec", 0);
558  vector< double > pool1Vec;
559  Field< double >::getVec( pool1, "n", pool1Vec );
560  assert( pool1Vec.size() == ndc );
561 
562  vector< double > pool2Vec;
563  Field< double >::getVec( pool2, "n", pool2Vec );
564  assert( pool2Vec.size() == ndc );
565 
566  vector< double > pool3Vec;
567  Field< double >::getVec( pool3, "n", pool3Vec );
568  assert( pool3Vec.size() == ndc );
569  double myTot1 = 0;
570  double myTot2 = 0;
571  double myTot3 = 0;
572  for ( unsigned int i = 0; i < nvec.size(); ++i ) {
573  assert( doubleEq( pool1Vec[i], nvec[i] ) );
574  myTot1 += nvec[i];
575  myTot2 += pool2Vec[i];
576  myTot3 += pool3Vec[i];
577  }
578  assert( doubleEq( myTot1, 1.0 ) );
579  assert( doubleEq( myTot2, 2.0 ) );
580  assert( doubleEq( myTot3, 3.0 ) );
581  assert( doubleEq( pool3Vec[0], 3.0 ) );
582 
583  /*
584  cout << "Small cell: " << endl;
585  for ( unsigned int i = 0; i < nvec.size(); ++i )
586  cout << nvec[i] << ", " << pool2Vec[i] << endl;
587  cout << endl;
588  */
589 
590 
591  s->doDelete( model );
592  cout << "." << flush;
593 }
594 
595 void testCellDiffn()
596 {
597  Id makeCompt( Id parentCompt, Id parentObj,
598  string name, double len, double dia, double theta );
599  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
600  double len = 40e-6;
601  double dia = 10e-6;
602  double diffLength = 1e-6;
603  double dt = 1.0e-1;
604  double runtime = 100.0;
605  double diffConst = 1.0e-12;
606  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
607  Id soma = makeCompt( Id(), model, "soma", dia, dia, 90 );
608  Id dend = makeCompt( soma, model, "dend", len, 3e-6, 0 );
609  Id branch1 = makeCompt( dend, model, "branch1", len, 2e-6, 45.0 );
610  Id branch2 = makeCompt( dend, model, "branch2", len, 2e-6, -45.0 );
611  Id twig1 = makeCompt( branch1, model, "twig1", len, 1.5e-6, 90.0 );
612  Id twig2 = makeCompt( branch1, model, "twig2", len, 1.5e-6, 0.0 );
613 
614  Id nm = s->doCreate( "NeuroMesh", model, "neuromesh", 1 );
615  Field< double >::set( nm, "diffLength", diffLength );
616  Field< string >::set( nm, "geometryPolicy", "cylinder" );
617  Field< string >::set( nm, "subTreePath", "/model/#" );
618  //SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );
619  unsigned int ns = Field< unsigned int >::get( nm, "numSegments" );
620  assert( ns == 6 );
621  unsigned int ndc = Field< unsigned int >::get( nm, "numDiffCompts" );
622  assert( ndc == 210 );
623  Id pool1 = s->doCreate( "Pool", nm, "pool1", 1 );
624  Field< double >::set( pool1, "diffConst", diffConst );
625  Id pool2 = s->doCreate( "Pool", nm, "pool2", 1 );
626  Field< double >::set( pool2, "diffConst", diffConst );
627 
628  Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
629  Field< Id >::set( dsolve, "compartment", nm );
630  s->doUseClock( "/model/dsolve", "process", 1 );
631  s->doSetClock( 1, dt );
632  // Next: build by setting path
633  Field< string >::set( dsolve, "path", "/model/neuromesh/pool#" );
634 
635  vector< double > nvec =
636  LookupField< unsigned int, vector< double > >::get(
637  dsolve, "nVec", 0);
638  assert( nvec.size() == ndc );
639  assert( pool1.element()->numData() == ndc );
640  Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
641  Field< double >::set( ObjId( pool2, ndc - 1 ), "nInit", 2.0 );
642 
643  s->doReinit();
644  s->doStart( runtime );
645 
646  nvec = LookupField< unsigned int, vector< double > >::get(
647  dsolve, "nVec", 0);
648  vector< double > pool1Vec;
649  Field< double >::getVec( pool1, "n", pool1Vec );
650  assert( pool1Vec.size() == ndc );
651 
652  vector< double > pool2Vec;
653  Field< double >::getVec( pool2, "n", pool2Vec );
654  assert( pool2Vec.size() == ndc );
655  double myTot1 = 0;
656  double myTot2 = 0;
657  for ( unsigned int i = 0; i < nvec.size(); ++i ) {
658  assert( doubleEq( pool1Vec[i], nvec[i] ) );
659  myTot1 += nvec[i];
660  myTot2 += pool2Vec[i];
661  }
662  assert( doubleEq( myTot1, 1.0 ) );
663  assert( doubleEq( myTot2, 2.0 ) );
664 
665  /*
666  cout << endl;
667  cout << "Big cell: " << endl;
668  for ( unsigned int i = 0; i < nvec.size(); ++i )
669  cout << nvec[i] << ", " << pool2Vec[i] << endl;
670  cout << endl;
671  */
672 
673 
674  s->doDelete( model );
675  cout << "." << flush;
676 }
677 
678 void testCylDiffnWithStoich()
679 {
680  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
681  double len = 25e-6;
682  double r0 = 1e-6;
683  double r1 = 1e-6;
684  double diffLength = 1e-6; // 1e-6 is the highest dx for which error is OK
685  double runtime = 10.0;
686  double dt0 = 0.1; // Used for diffusion. 0.2 is the highest dt for which the error is in bounds
687  double dt1 = 1; // Used for chem.
688  double diffConst = 1.0e-12;
689  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
690  Id cyl = s->doCreate( "CylMesh", model, "cyl", 1 );
691  Field< double >::set( cyl, "r0", r0 );
692  Field< double >::set( cyl, "r1", r1 );
693  Field< double >::set( cyl, "x0", 0 );
694  Field< double >::set( cyl, "x1", len );
695  Field< double >::set( cyl, "diffLength", diffLength );
696  unsigned int ndc = Field< unsigned int >::get( cyl, "numMesh" );
697  assert( ndc == static_cast< unsigned int >( round( len / diffLength )));
698  Id pool1 = s->doCreate( "Pool", cyl, "pool1", 1 );
699  Id pool2 = s->doCreate( "Pool", cyl, "pool2", 1 );
700  Field< double >::set( pool1, "diffConst", diffConst );
701  Field< double >::set( pool2, "diffConst", diffConst/2 );
702 
703  Id stoich = s->doCreate( "Stoich", model, "stoich", 1 );
704  Id ksolve = s->doCreate( "Ksolve", model, "ksolve", 1 );
705  Id dsolve = s->doCreate( "Dsolve", model, "dsolve", 1 );
706  Field< Id >::set( stoich, "compartment", cyl );
707  Field< Id >::set( stoich, "ksolve", ksolve );
708  Field< Id >::set( stoich, "dsolve", dsolve );
709  Field< string >::set( stoich, "path", "/model/cyl/#" );
710  assert( pool1.element()->numData() == ndc );
711 
712  // Then find a way to test it.
713  vector< double > poolVec;
714  Field< double >::set( ObjId( pool1, 0 ), "nInit", 1.0 );
715  Field< double >::set( ObjId( pool2, 0 ), "nInit", 1.0 );
716  Field< double >::getVec( pool1, "nInit", poolVec );
717  assert( poolVec.size() == ndc );
718  assert( doubleEq( poolVec[0], 1.0 ) );
719  assert( doubleEq( poolVec[1], 0.0 ) );
720 
721  vector< double > nvec =
722  LookupField< unsigned int, vector< double > >::get(
723  dsolve, "nVec", 0);
724  assert( nvec.size() == ndc );
725 
726  // Next: build by doing reinit
727  s->doUseClock( "/model/dsolve", "process", 0 );
728  s->doUseClock( "/model/ksolve", "process", 1 );
729  s->doSetClock( 0, dt0 );
730  s->doSetClock( 1, dt1 );
731  s->doReinit();
732  s->doStart( runtime );
733 
734  nvec = LookupField< unsigned int, vector< double > >::get(
735  dsolve, "nVec", 0);
736  Field< double >::getVec( pool1, "n", poolVec );
737  assert( nvec.size() == poolVec.size() );
738  for ( unsigned int i = 0; i < nvec.size(); ++i )
739  assert( doubleEq( nvec[i], poolVec[i] ) );
740  /*
741  cout << endl;
742  for ( unsigned int i = 0; i < nvec.size(); ++i )
743  cout << nvec[i] << " ";
744  cout << endl;
745  */
746 
747  double dx = diffLength;
748  double err = 0.0;
749  double analyticTot = 0.0;
750  double myTot = 0.0;
751  for ( unsigned int i = 0; i < nvec.size(); ++i ) {
752  double x = i * dx + dx * 0.5;
753  // This part is the solution as a func of x,t.
754  double y = dx * // This part represents the init n of 1 in dx
755  ( 1.0 / sqrt( PI * diffConst * runtime ) ) *
756  exp( -x * x / ( 4 * diffConst * runtime ) );
757  err += ( y - nvec[i] ) * ( y - nvec[i] );
758  //cout << i << " " << x << " " << y << " " << conc[i] << endl;
759  analyticTot += y;
760  myTot += nvec[i];
761  }
762  assert( doubleEq( myTot, 1.0 ) );
763  // cout << "analyticTot= " << analyticTot << ", myTot= " << myTot << endl;
764  assert( err < 1.0e-5 );
765 
766 
767  s->doDelete( model );
768  cout << "." << flush;
769 }
770 #if 0
771 void testBuildTree()
772 {
773  static double test[] = {
774  1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0,
775  5, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0,
776  7, 0, 8, 9, 0, 0, 0, 0, 10, 0, 0,
777  11, 0, 12, 13, 0, 0, 0, 0, 0, 0, 14,
778  0, 0, 0, 0, 15, 16, 0, 0, 17, 18, 0,
779  0, 0, 0, 0, 19, 20, 21, 0, 22, 0, 0,
780  0, 0, 0, 0, 0, 23, 24, 25, 0, 0, 0,
781  0, 0, 0, 0, 0, 0, 26, 27, 0, 0, 0,
782  0, 0, 28, 0, 29, 30, 0, 0, 31, 0, 0,
783  0, 0, 0, 0, 32, 0, 0, 0, 0, 33, 0,
784  0, 0, 0, 34, 0, 0, 0, 0, 0, 0, 35,
785  };
786  const unsigned int numCompts = 11;
787  FastMatrixElim fe;
788  fe.makeTestMatrix( test, numCompts );
789  fe.print();
790  vector< unsigned int > parentVoxel;
791  bool ret = fe.buildTree( 0, parentVoxel );
792  assert( ret );
793  assert( parentVoxel[0] == static_cast< unsigned int >( -1 ) );
794  assert( parentVoxel[1] == 0 );
795  assert( parentVoxel[2] == 0 );
796  assert( parentVoxel[3] == 0 );
797  assert( parentVoxel[8] == 2 );
798  assert( parentVoxel[10] == 3 );
799  assert( parentVoxel[5] == 4 );
800  assert( parentVoxel[9] == 4 );
801  assert( parentVoxel[6] == 5 );
802  assert( parentVoxel[7] == 6 );
803 
804  assert( parentVoxel[10] == 2 );
805 
806  /*
807  * This is the sequence of traversal. x means empty, s means sibling,
808  * . means below diagonal. The numbers are the sequence.
809  static double traverseIndex[] = {
810  // col 1 2 3 4 5 6 7 8 9 10
811  #, 1, 2, 3, x, x, x, x, x, x, x,
812  ., #, x, x, x, x, x, x, x, x, x,
813  ., x, #, s, x, x, x, x, 4, x, x,
814  ., x, ., #, x, x, x, x, x, x, 5,
815  x, x, x, x, #, 6, x, x, s, 7, x,
816  x, x, x, x, ., #, 8, x, s, x, x,
817  x, x, x, x, x, ., #, 9, x, x, x,
818  x, x, x, x, x, x, ., #, x, x, x,
819  x, x, ., x, ., ., x, x, #, x, x,
820  x, x, x, x, ., x, x, x, x, #, x,
821  x, x, x, ., x, x, x, x, x, x, #,
822  };
823  */
824 }
825 #endif
826 
827 void testCalcJunction()
828 {
829  Shell* s = reinterpret_cast< Shell* >( Id().eref().data() );
830  // Make a neuron with same-size dend and spine. PSD is tiny.
831  // Put a, b, c in dend, b, c, d in spine, c, d, f in psd. No reacs.
832  // See settling of all concs by diffusion, pairwise.
833  Id model = s->doCreate( "Neutral", Id(), "model", 1 );
834  Id dend = s->doCreate( "Compartment", model, "dend", 1 );
835  Id neck = s->doCreate( "Compartment", model, "spine_neck", 1 );
836  Id head = s->doCreate( "Compartment", model, "spine_head", 1 );
837  Field< double >::set( dend, "x", 10e-6 );
838  Field< double >::set( dend, "diameter", 2e-6 );
839  Field< double >::set( dend, "length", 10e-6 );
840  Field< double >::set( neck, "x0", 9e-6 );
841  Field< double >::set( neck, "x", 9e-6 );
842  Field< double >::set( neck, "y", 1e-6 );
843  Field< double >::set( neck, "diameter", 0.5e-6 );
844  Field< double >::set( neck, "length", 1.0e-6 );
845  Field< double >::set( head, "x0", 9e-6 );
846  Field< double >::set( head, "x", 9e-6 );
847  Field< double >::set( head, "y0", 1e-6 );
848  Field< double >::set( head, "y", 11e-6 );
849  Field< double >::set( head, "diameter", 2e-6 );
850  Field< double >::set( head, "length", 10e-6 );
851  s->doAddMsg( "Single", ObjId( dend ), "raxial", ObjId( neck ), "axial");
852  s->doAddMsg( "Single", ObjId( neck ), "raxial", ObjId( head ), "axial");
853 
854  Id nm = s->doCreate( "NeuroMesh", model, "nm", 1 );
855  Field< double >::set( nm, "diffLength", 10e-6 );
856  Field< bool >::set( nm, "separateSpines", true );
857  Id sm = s->doCreate( "SpineMesh", model, "sm", 1 );
858  Id pm = s->doCreate( "PsdMesh", model, "pm", 1 );
859  ObjId mid = s->doAddMsg( "Single", ObjId( nm ), "spineListOut", ObjId( sm ), "spineList" );
860  assert( !mid.bad() );
861  mid = s->doAddMsg( "Single", ObjId( nm ), "psdListOut", ObjId( pm ), "psdList" );
862  Field< string >::set( nm, "subTreePath", "/model/#" );
863  //SetGet2< Id, string >::set( nm, "cellPortion", model, "/model/#" );
864 
865  vector< Id > pools( 9 );
866  static string names[] = {"a", "b", "c", "b", "c", "d", "c", "d", "e" };
867  static Id parents[] = {nm, nm, nm, sm, sm, sm, pm, pm, pm};
868  for ( unsigned int i = 0; i < 9; ++i ) {
869  pools[i] = s->doCreate( "Pool", parents[i], names[i], 1 );
870  assert( pools[i] != Id() );
871  Field< double >::set( pools[i], "concInit", 1.0 + 1.0 * i );
872  Field< double >::set( pools[i], "diffConst", 1e-11 );
873  if ( i < 6 ) {
874  double vol = Field< double >::get( pools[i], "volume" );
875  assert( doubleEq( vol, 10e-6 * 1e-12 * PI ) );
876  }
877  }
878  Id dendsolve = s->doCreate( "Dsolve", model, "dendsolve", 1 );
879  Id spinesolve = s->doCreate( "Dsolve", model, "spinesolve", 1 );
880  Id psdsolve = s->doCreate( "Dsolve", model, "psdsolve", 1 );
881  Field< Id >::set( dendsolve, "compartment", nm );
882  Field< Id >::set( spinesolve, "compartment", sm );
883  Field< Id >::set( psdsolve, "compartment", pm );
884  Field< string >::set( dendsolve, "path", "/model/nm/#" );
885  Field< string >::set( spinesolve, "path", "/model/sm/#" );
886  Field< string >::set( psdsolve, "path", "/model/pm/#" );
887  assert( Field< unsigned int >::get( dendsolve, "numAllVoxels" ) == 1 );
888  assert( Field< unsigned int >::get( spinesolve, "numAllVoxels" ) == 1 );
889  assert( Field< unsigned int >::get( psdsolve, "numAllVoxels" ) == 1 );
890  assert( Field< unsigned int >::get( dendsolve, "numPools" ) == 3 );
891  assert( Field< unsigned int >::get( spinesolve, "numPools" ) == 3 );
892  assert( Field< unsigned int >::get( psdsolve, "numPools" ) == 3 );
893  SetGet2< Id, Id >::set( dendsolve, "buildNeuroMeshJunctions",
894  spinesolve, psdsolve );
895  s->doSetClock( 0, 0.01 );
896  s->doUseClock( "/model/#solve", "process", 0 );
897  s->doReinit();
898  s->doStart( 100 );
899 
900  /*
901  for ( unsigned int i = 0; i < 9; ++i ) {
902  double c = Field< double >::get( pools[i], "conc" );
903  double n = Field< double >::get( pools[i], "n" );
904  double v = Field< double >::get( pools[i], "volume" );
905  cout << pools[i].path() << ": " << c << ", " << n << ", " <<
906  n / v << ", " <<
907  v << endl;
908  }
909  */
910  s->doDelete( model );
911  cout << "." << flush;
912 }
913 
914 void testDiffusion()
915 {
916  testSorting();
917  testFastMatrixElim();
918  testSetDiffusionAndTransport();
919  testCylDiffn();
920  testTaperingCylDiffn();
921  testSmallCellDiffn();
922  testCellDiffn();
923  testCylDiffnWithStoich();
924  testCalcJunction();
925 }
Shell::doStart
void doStart(double runtime, bool notify=false)
Definition: Shell.cpp:332
Eref::data
char * data() const
Definition: Eref.cpp:41
SparseMatrix::get
T get(unsigned int row, unsigned int column) const
Definition: SparseMatrix.h:253
makeCompt
Id makeCompt(Id parentCompt, Id parentObj, string name, double len, double dia, double theta)
Definition: testMesh.cpp:855
Id::element
Element * element() const
Synonym for Id::operator()()
Definition: Id.cpp:113
FastMatrixElim::buildBackwardSub
void buildBackwardSub(vector< unsigned int > &diag, vector< Triplet< double > > &bops, vector< double > &diagVal)
Definition: FastMatrixElim.cpp:342
FastMatrixElim::advance
static void advance(vector< double > &y, const vector< Triplet< double > > &ops, const vector< double > &diagVal)
Definition: FastMatrixElim.cpp:472
Shell::doSetClock
void doSetClock(unsigned int tickNum, double dt)
Definition: Shell.cpp:377
ObjId::bad
bool bad() const
Definition: ObjId.cpp:18
Field
Definition: SetGet.h:236
testDiffusion
void testDiffusion()
Definition: testDiffusion.cpp:914
testCylDiffnWithStoich
void testCylDiffnWithStoich()
Definition: testDiffusion.cpp:678
ObjId
Definition: ObjId.h:20
FastMatrixElim::hinesReorder
bool hinesReorder(const vector< unsigned int > &parentVoxel, vector< unsigned int > &lookupOldRowsFromNew)
Definition: FastMatrixElim.cpp:139
Id::eref
Eref eref() const
Definition: Id.cpp:125
Field::set
static bool set(const ObjId &dest, const string &field, A arg)
Definition: SetGet.h:245
Shell::doCreate
Id doCreate(string type, ObjId parent, string name, unsigned int numData, NodePolicy nodePolicy=MooseBlockBalance, unsigned int preferredNode=1)
Definition: Shell.cpp:181
testTaperingCylDiffn
void testTaperingCylDiffn()
Definition: testDiffusion.cpp:422
testSmallCellDiffn
void testSmallCellDiffn()
Definition: testDiffusion.cpp:498
testFastMatrixElim
void testFastMatrixElim()
Definition: testDiffusion.cpp:46
doubleEq
bool doubleEq(double x, double y)
Definition: doubleEq.cpp:16
testCellDiffn
void testCellDiffn()
Definition: testDiffusion.cpp:595
testCylDiffn
void testCylDiffn()
Definition: testDiffusion.cpp:341
Element::numData
virtual unsigned int numData() const =0
Returns number of data entries across all nodes.
testSetDiffusionAndTransport
void testSetDiffusionAndTransport()
Definition: testDiffusion.cpp:288
SparseMatrix::print
void print() const
Definition: SparseMatrix.h:650
FastMatrixElim::makeTestMatrix
void makeTestMatrix(const double *test, unsigned int numCompts)
Definition: FastMatrixElim.cpp:244
Shell::doReinit
void doReinit()
Definition: Shell.cpp:362
name
static char name[]
Definition: mfield.cpp:401
Shell::doAddMsg
ObjId doAddMsg(const string &msgType, ObjId src, const string &srcField, ObjId dest, const string &destField)
Definition: Shell.cpp:269
FastMatrixElim::setDiffusionAndTransport
void setDiffusionAndTransport(const vector< unsigned int > &parentVoxel, double diffConst, double motorConst, double dt)
Definition: FastMatrixElim.cpp:405
testSorting
void testSorting()
Definition: testDiffusion.cpp:253
PI
const double PI
Definition: consts.cpp:12
Shell::doDelete
bool doDelete(ObjId oid)
Definition: Shell.cpp:259
Id
Definition: Id.h:17
checkAns
double checkAns(const double *m, unsigned int numCompts, const double *ans, const double *rhs)
Definition: testDiffusion.cpp:29
Shell::doUseClock
void doUseClock(string path, string field, unsigned int tick)
Definition: Shell.cpp:382
Field::get
static A get(const ObjId &dest, const string &field)
Definition: SetGet.h:284
Field::getVec
static void getVec(ObjId dest, const string &field, vector< A > &vec)
Definition: SetGet.h:317
testCalcJunction
void testCalcJunction()
Definition: testDiffusion.cpp:827
FastMatrixElim::buildForwardElim
void buildForwardElim(vector< unsigned int > &diag, vector< Triplet< double > > &fops)
Definition: FastMatrixElim.cpp:275
Shell
Definition: Shell.h:43
sortByColumn
void sortByColumn(vector< unsigned int > &col, vector< double > &entry)
Definition: FastMatrixElim.cpp:223
SetGet2::set
static bool set(const ObjId &dest, const string &field, A1 arg1, A2 arg2)
Definition: SetGet.h:365