UltraScan III
us_lamm_astfvm.cpp
Go to the documentation of this file.
1 
3 #include "us_lamm_astfvm.h"
4 #include "us_math2.h"
5 #include "us_constants.h"
6 #include "us_astfem_rsa.h"
7 #include "us_settings.h"
8 #include "us_dataIO.h"
9 
11 //
12 // Mesh
13 //
15 US_LammAstfvm::Mesh::Mesh( double xl, double xr, int Nelem, int Opt )
16 {
17  int i;
18  dbg_level = US_Settings::us_debug();
19 
20  // constants
21  MaxRefLev = 20;
22  MonScale = 1;
23  MonCutoff = 1000;
24  SmoothingWt = 0.7;
25  SmoothingCyl = 4;
26 
27  Ne = Nelem;
28  Nv = Ne + 1;
29  Eid = new int [ Ne ];
30  RefLev = new int [ Ne ];
31  Mark = new int [ Ne ];
32  MeshDen = new double [ Ne ];
33  x = new double [ Nv ];
34 
35  // uniform
36  if ( Opt == 0 )
37  {
38  for ( i = 0; i < Nv; i++ )
39  {
40  x[ i ] = xl + (double)i / (double)Ne * ( xr - xl );
41  }
42 
43  for ( i = 0; i < Ne; i += 2 ) Eid[ i ] = 1 ;
44  for ( i = 1; i < Ne; i += 2 ) Eid[ i ] = 4 ;
45  for ( i = 0; i < Ne; i++ ) RefLev[ i ] = 0;
46  }
47 }
48 
50 //
51 // ~Mesh
52 //
54 US_LammAstfvm::Mesh::~Mesh()
55 {
56  delete [] x;
57  delete [] Eid;
58  delete [] RefLev;
59  delete [] MeshDen;
60  delete [] Mark;
61 }
62 
64 //
65 // ComputeMeshDen_D3
66 //
68 void US_LammAstfvm::Mesh::ComputeMeshDen_D3( double *u0, double *u1 )
69 {
70  int i;
71  int i2;
72  double h;
73 
74  double* D20 = new double [ Ne ];
75  double* D21 = new double [ Ne ];
76  double* D30 = new double [ Nv ];
77  double* D31 = new double [ Nv ];
78  // 2nd derivative on elems
79  for ( i = 0; i < Ne; i++ )
80  {
81  h = ( x[ i + 1 ] - x[ i ] ) / 2;
82  i2 = i * 2;
83  D20[ i ] = ( u0[ i2 + 2 ] - 2 * u0[ i2 + 1 ] + u0[ i2 ] ) / ( h * h );
84  D21[ i ] = ( u1[ i2 + 2 ] - 2 * u1[ i2 + 1 ] + u1[ i2 ] ) / ( h * h );
85  }
86 
87  // 3rd derivative at nodes
88  for ( i = 1; i < Nv - 1; i++ )
89  {
90  h = ( x[ i + 1 ] - x[ i - 1 ] ) / 2;
91  D30[ i ] = ( D20[ i ] - D20[ i - 1 ] ) / h;
92  D31[ i ] = ( D21[ i ] - D21[ i - 1 ] ) / h;
93  }
94 
95  // 3rd derivative on elems = average D3u at nodes
96  // here use D2 to store 3rd order derivatives
97  for ( i = 1; i < Ne - 1; i++ )
98  {
99  D20[ i ] = ( D30[ i + 1 ] + D30[ i ] ) / 2;
100  D21[ i ] = ( D31[ i + 1 ] + D31[ i ] ) / 2;
101  }
102 
103  D20[ 0 ] = D20[ 1 ]; // Extropolate 3rd derivative to end pts
104  D21[ 0 ] = D21[ 1 ];
105  D20[ Ne - 1 ] = D20[ Ne - 2 ];
106  D21[ Ne - 1 ] = D21[ Ne - 2 ];
107 
108  // level-off and smoothing
109  for ( i = 0; i < Ne; i++ )
110  {
111  MeshDen[ i ] = pow( 1 + MonScale
112  * ( fabs( D20[ i ] ) + fabs( D21[ i ] ) ), 0.33333 );
113 
114  if ( MeshDen[ i ] > MonCutoff ) MeshDen[ i ] = MonCutoff;
115  }
116 
117  Smoothing( Ne, MeshDen, SmoothingWt, SmoothingCyl );
118 
119  delete [] D20;
120  delete [] D21;
121  delete [] D30;
122  delete [] D31;
123 }
124 
125 
126 
127 
129 // y[i] = (1-Wt)/2 * y1 +
130 // Smoothing: Wt * y2 + // for y1,y2,y3 unsmoothed points
131 // (1-Wt)/2 * y3; // around y[i]
133 void US_LammAstfvm::Mesh::Smoothing( int n, double *y, double Wt, int Cycle )
134 {
135  int s;
136  int i;
137  double Wt1;
138  double Wt2;
139  Wt1 = 1. - Wt; // sum of outside pt. weights
140  Wt2 = Wt1 * 0.5; // weight, each outside pt.
141 
142  for ( s = 0; s < Cycle; s++ )
143  {
144  double y1;
145  double y2 = y[ 0 ];
146  double y3 = y[ 1 ];
147 
148  y[ 0 ] = Wt * y3 + Wt1 * y2; // 1st smoothed point
149 
150  for ( i = 1; i < n - 1; i++ )
151  { // smooth all points except first and last
152  y1 = y2; // cycle the 3 points
153  y2 = y3; // around the current one
154  y3 = y[ i + 1 ];
155 
156  y[ i ] = Wt * y2 + Wt2 * ( y1 + y3 ); // smoothed point
157  }
158 
159  y[ n - 1 ] = Wt * y2 + Wt1 * y3; // last smoothed point
160  }
161 }
162 
164 //
165 // Unrefine
166 //
168 void US_LammAstfvm::Mesh::Unrefine( double alpha )
169 {
170  int i;
171  int i1;
172  int Nv1;
173  int Ne1;
174  double* x1;
175  int* Eid1; // elem id
176  int* RefLev1; // refinement level of an elem
177  double* MeshDen1; // desired mesh density
178 
179  while( 1 )
180  {
181  // set unref mark on each elem
182  for ( i = 0; i < Ne; i++ ) Mark[ i ] = 0;
183 
184  Ne1 = Ne;
185 
186  for ( i = 0; i < Ne - 1; i++ )
187  {
188  if ( RefLev[ i ] == RefLev[ i + 1 ] &&
189  ( Eid[ i ] / 2 ) == ( Eid[ i + 1 ] / 2 ) &&
190  ( x[ i + 1 ] - x[ i ] ) * MeshDen[ i ] < alpha &&
191  ( x[ i + 2 ] - x[ i + 1 ] ) * MeshDen[ i + 1 ] < alpha )
192  {
193  Mark[ i ] = 1;
194  Mark[ i + 1 ] = 1;
195  Ne1 --;
196  i ++;
197  }
198  }
199 
200  if ( Ne1 == Ne ) return; // no more unrefine
201 
202  // reallocate memory for the new mesh
203  Nv1 = Ne1 + 1;
204  Eid1 = new int [ Ne1 ];
205  RefLev1 = new int [ Ne1 ];
206  MeshDen1 = new double [ Ne1 ];
207  x1 = new double [ Nv1 ];
208 
209  // loop to combine eligible elem pairs
210  x1[ 0 ] = x[ 0 ];
211  i1 = 0;
212 
213  for ( i = 0; i < Ne; i++ )
214  {
215  if ( Mark[ i ] == 1 && Mark[ i + 1 ] == 1 )
216  { // combine two elems
217  x1[ i1 + 1 ] = x[ i + 2 ];
218  Eid1[ i1 ] = Eid[ i ] / 2;
219  RefLev1[ i1 ] = RefLev[ i ] - 1;
220  MeshDen1[ i1 ] = ( MeshDen[ i ] + MeshDen[ i + 1 ] ) / 2;
221  i1++;
222  i++;
223  }
224 
225  else
226  { // no change
227  x1[ i1 + 1 ] = x[ i + 1 ];
228  Eid1[ i1 ] = Eid[ i ];
229  RefLev1[ i1 ] = RefLev[ i ];
230  MeshDen1[ i1 ] = MeshDen[ i ];
231  i1++;
232  }
233  }
234 
235  // delete memory for old mesh
236  delete [] x;
237  delete [] Eid;
238  delete [] RefLev;
239  delete [] MeshDen;
240  delete [] Mark;
241 
242  Ne = Ne1;
243  Nv = Nv1;
244  x = x1;
245  Eid = Eid1;
246  MeshDen = MeshDen1;
247  RefLev = RefLev1;
248  Mark = new int [ Ne1 ];
249 
250  } // while
251 
252 }
253 
254 
256 //
257 // Refine
258 //
260 void US_LammAstfvm::Mesh::Refine( double beta )
261 {
262  int k;
263  int ke;
264  int Nv1;
265  int Ne1; // number of elements
266  int* Eid1; // element id
267  int* RefLev1; // refinement level of an element
268  double* MeshDen1; // desired mesh density
269  double* x1;
270 
271  while( 1 )
272  {
273 
274  // set marker for elements that need to be refined
275  for ( k = 0; k < Ne; k++ )
276  {
277  if ( ( x[ k + 1 ] - x[ k ] ) * MeshDen[ k ] > beta &&
278  RefLev[ k ] < MaxRefLev )
279  Mark[ k ] = 1;
280 
281  else
282  Mark[ k ] = 0;
283  }
284 
285  for ( k = 0; k < Ne - 1; k++ )
286  { // RefLev differs at most 2 for nabos
287  int rldiff = RefLev[ k ] - RefLev[ k + 1 ];
288 
289  if ( rldiff < (-1) )
290  Mark[ k ] = 1;
291 
292  else if ( rldiff > 1 )
293  Mark[ k + 1 ] = 1;
294  }
295 
296  Ne1 = Ne;
297 
298  for ( k = 0; k < Ne; k++ )
299  if ( Mark[ k ] == 1 ) Ne1++;
300 
301  if ( Ne1 == Ne ) return; // no more elements need refine
302 
303  // allocate memory for new mesh
304  Nv1 = Ne1 + 1;
305  x1 = new double [ Nv1 ];
306  Eid1 = new int [ Ne1 ];
307  RefLev1 = new int [ Ne1 ];
308  MeshDen1 = new double [ Ne1 ];
309 
310  ke = 0;
311  x1[ 0 ] = x[ 0 ];
312 
313  for ( k = 0; k < Ne; k++ )
314  {
315  if ( Mark[ k ] == 0 )
316  { // no refine on elem-k
317  x1[ ke + 1 ] = x[ k + 1 ];
318  Eid1[ ke ] = Eid[ k ];
319  RefLev1[ ke ] = RefLev[ k ];
320  MeshDen1[ ke ] = MeshDen[ k ];
321  ke++;
322  }
323 
324  else
325  { // refine k-th elem
326  x1[ ke + 1 ] = ( x[ k ] + x[ k + 1 ] ) / 2;
327  Eid1[ ke ] = Eid[ k ] * 2;
328  RefLev1[ ke ] = RefLev[ k ] + 1;
329  MeshDen1[ ke ] = MeshDen[ k ];
330 
331  x1[ ke + 2 ] = x[ k + 1 ];
332  Eid1[ ke + 1 ] = Eid1[ ke ] + 1;
333  RefLev1[ ke + 1 ] = RefLev1[ ke ];
334  MeshDen1[ ke + 1 ] = MeshDen[ k ];
335  ke += 2;
336  } // if
337  } // for
338 
339  delete [] x;
340  delete [] Eid;
341  delete [] RefLev;
342  delete [] Mark;
343  delete [] MeshDen;
344 
345  Ne = Ne1;
346  Nv = Nv1;
347  x = x1;
348  Eid = Eid1;
349  RefLev = RefLev1;
350  MeshDen = MeshDen1;
351  Mark = new int [ Ne1 ];
352  } // while
353 }
354 
356 //
357 // RefineMesh
358 //
360 void US_LammAstfvm::Mesh::RefineMesh( double *u0, double *u1, double ErrTol )
361 {
362  const double sqrt3 = sqrt( 3.0 );
363  const double onethird = 1.0 / 3.0;
364  // refinement threshhold: h*|D_3u|^(1/3) > beta
365  //double beta = pow( ErrTol * 6 / ( 2 * sqrt( 3 ) / 72 ), 1. / 3. );
366  // Simplify the above:
367  double beta = 6.0 * pow( ErrTol / sqrt3, onethird );
368 
369  // coarsening threshhold: h*|D_3u|^(1/3) < alpha
370  double alpha = beta / 4;
371 //DbgLv(3) << "RefMesh: beta alpha" << beta << alpha;
372 
373  ComputeMeshDen_D3( u0, u1 );
374  Smoothing ( Ne, MeshDen, 0.7, 4 );
375  Unrefine ( alpha );
376  Refine ( beta );
377 }
378 
380 //
381 // InitMesh
382 //
384 void US_LammAstfvm::Mesh::InitMesh( double s, double D, double w2 )
385 {
386  int j;
387  double D0;
388  double nu0;
389  double nu1;
390  double nu;
391  double t;
392  double m2;
393  double b2;
394  double x2;
395  double* u0;
396  double* u1;
397 
398  D0 = 1.e-4;
399  nu0 = s * w2 / D0;
400  nu1 = s * w2 / D;
401 
402  m2 = x[ 0 ] * x[ 0 ];
403  b2 = x[ Ne ] * x[ Ne ];
404 
405  // FILE *fout;
406  // fout = fopen("ti.tmp", "w");
407 
408  for ( t = 0; t < 1; t = t + 0.1 )
409  {
410  u0 = new double [ 2 * Nv - 1 ];
411  u1 = new double [ 2 * Nv - 1 ];
412 
413  nu = pow( nu1, t ) * pow( nu0, 1 - t );
414 
415  for ( j = 0; j < Nv; j++ )
416  {
417  x2 = x[ j ] * x[ j ];
418  u0[ 2 * j ] = exp( nu * ( m2 - x2 ) ) * ( nu * nu * m2 );
419  u1[ 2 * j ] = exp( nu * ( x2 - b2 ) ) * ( nu * nu * b2 );
420  }
421 
422  for ( j = 0; j < Ne; j++ )
423  {
424  x2 = ( x[ j ] + x[ j + 1 ] ) / 2;
425  x2 = x2 * x2;
426  u0[ 2 * j + 1 ] = exp( nu * ( m2 - x2 ) ) * ( nu * nu * m2 );
427  u1[ 2 * j + 1 ] = exp( nu * ( x2 - b2 ) ) * ( nu * nu * b2 );
428  }
429 
430  RefineMesh( u0, u1, 1.e-4 );
431 
432  delete [] u0;
433  delete [] u1;
434  }
435 }
436 
437 // create salt data set by solving ideal astfem equations
438 US_LammAstfvm::SaltData::SaltData( US_Model amodel,
439  US_SimulationParameters asparms,
440  US_DataIO::RawData* asim_data )
441 {
442  dbg_level = US_Settings::us_debug();
443  model = amodel;
444  simparms = asparms;
445  sa_data = *asim_data;
446 DbgLv(2) << "SaltD: sa_data avg.temp." << sa_data.average_temperature();
447 DbgLv(2) << "SaltD: asim_data avg.temp." << asim_data->average_temperature();
448 
449  Nt = sa_data.scanCount();
450  Nx = sa_data.pointCount();
451 
452  model.components.resize( 1 );
453  model.components[ 0 ] = amodel.components[ amodel.coSedSolute ];
454  model.coSedSolute = -1;
455  // use salt's molar concentration, if possible, for initial concentration
456  double conc0 = model.components[ 0 ].molar_concentration;
457  conc0 = ( conc0 < 0.01 ) ?
458  model.components[ 0 ].signal_concentration : conc0;
459  conc0 = ( conc0 < 0.01 ) ? 2.5 : conc0;
460 //conc0=3100.0;
461  model.components[ 0 ].signal_concentration = conc0;
462 
463  simparms.meshType = US_SimulationParameters::ASTFEM;
464  simparms.gridType = US_SimulationParameters::MOVING;
465 
466  simparms.radial_resolution =
467  ( sa_data.radius( Nx - 1 ) - sa_data.radius( 0 ) ) / (double)( Nx - 1 );
468  simparms.firstScanIsConcentration = false;
469 DbgLv(2) << "SaltD: Nx" << Nx << "r0 rn ri" << sa_data.radius( 0 )
470  << sa_data.radius( Nx - 1 ) << simparms.radial_resolution;
471 
472  US_Astfem_RSA* astfem = new US_Astfem_RSA( model, simparms );
473 
474  //astfem->setTimeInterpolation( true );
475  astfem->setTimeInterpolation( false );
476  astfem->set_debug_flag( dbg_level );
477 
478  for ( int ii = 0; ii < Nt; ii++ )
479  for ( int jj = 0; jj < Nx; jj++ )
480  sa_data.setValue( ii, jj, 0.0 );
481 
482 DbgLv(2) << "SaltD: model comps" << model.components.size();
483 DbgLv(2) << "SaltD: amodel comps" << amodel.components.size();
484 DbgLv(2) << "SaltD: comp0 s d s_conc" << model.components[0].s
485  << model.components[0].D << model.components[0].signal_concentration;
486 DbgLv(2) << "SaltD:fem: m b s D rpm" << simparms.meniscus << simparms.bottom
487  << model.components[0].s << model.components[0].D
488  << simparms.speed_step[0].rotorspeed;
489 DbgLv(2) << "SaltD: (0)Nt Nx" << Nt << Nx << "temperature" << sa_data.scanData[0].temperature;
490 if(dbg_level>0) {
491  qDebug() << "SaltD: model";
492  model.debug();
493  qDebug() << "SaltD: simparms:";
494  simparms.debug();
495 }
496 
497  astfem->set_simout_flag( true ); // set flag to output raw simulation
498 DbgLv(2) << "SaltD: (2)Nx" << Nx << "r0 r1 rm rn" << sa_data.radius(0)
499  << sa_data.radius(1) << sa_data.radius(Nx-2) << sa_data.radius(Nx-1);
500 
501  astfem->calculate( sa_data ); // solve equations to create data
502 
503  Nt = sa_data.scanCount();
504  Nx = sa_data.pointCount();
505 
506 DbgLv(2) << "SaltD: (3)Nx" << Nx << "r0 r1 rm rn" << sa_data.radius(0)
507  << sa_data.radius(1) << sa_data.radius(Nx-2) << sa_data.radius(Nx-1);
508 DbgLv(2) << "SaltD: Nt Nx" << Nt << Nx << "temp" << sa_data.scanData[0].temperature;
509 DbgLv(2) << "SaltD: sa sc0 sec omg" << sa_data.scanData[0].seconds
510  << sa_data.scanData[0].omega2t;
511  // Limit salt amplitudes to reasonable values
512  const double maxsalt = 1e100;
513  const double minsalt = -9e99;
514  const double minsala = 1e-100;
515  const double minsaln = -1e-100;
516  int nchg = 0;
517  for ( int ii = 0; ii < Nt; ii++ )
518  {
519  for ( int jj = 0; jj < Nx; jj++ )
520  {
521  double saltv = sa_data.value( ii, jj );
522  double salta = qAbs( saltv );
523 if(ii<2 && (jj+3)>Nx)
524 DbgLv(2) << "SaltD: ii jj" << ii << jj << "saltv salta" << saltv << salta;
525 
526  if ( saltv != saltv )
527  { // Nan!!!
528  nchg++;
529  if ( jj > 0 )
530  {
531  saltv = sa_data.value( ii, jj - 1 );
532  }
533  else
534  {
535  saltv = saltv > 0 ? maxsalt : minsalt;
536  }
537  sa_data.setValue( ii, jj, saltv );
538 DbgLv(2) << "SaltD: salt *NaN* ii jj adj-saltv" << ii << jj << saltv;
539  }
540 
541  else if ( salta > maxsalt || saltv < minsalt )
542  { // Amplitude too large
543  nchg++;
544  salta = qMin( maxsalt, qMax( minsalt, salta ) );
545  saltv = saltv > 0 ? salta : -salta;
546  sa_data.setValue( ii, jj, saltv );
547  }
548 
549  else if ( salta < minsala )
550  { // Amplitude too small
551  nchg++;
552  saltv = saltv > 0 ? minsala : minsaln;
553  sa_data.setValue( ii, jj, saltv );
554  }
555  }
556  }
557 DbgLv(2) << "SaltD: salt ampl limit changes" << nchg;
558 
559  //xs = new double [ Nx ];
560  //Cs0 = new double [ Nx ];
561  //Cs1 = new double [ Nx ];
562  xsVec .fill( 0.0, Nx );
563  Cs0Vec.fill( 0.0, Nx );
564  Cs1Vec.fill( 0.0, Nx );
565  xs = xsVec .data();
566  Cs0 = Cs0Vec.data();
567  Cs1 = Cs1Vec.data();
568 
569  if ( dbg_level > 2 )
570  { // save a copy of the salt data set so that it may be plotted for QC
571  QString safile = US_Settings::resultDir() + "/salt_data";
572  QDir dir;
573 
574  dir.mkpath( safile );
575  safile = safile + "/salt_data.RA.1.S.260.auc";
576  US_DataIO::writeRawData( safile, sa_data );
577  }
578 
579  delete astfem; // astfem solver no longer needed
580 
581  for ( int jj = 0; jj < Nx; jj++ )
582  { // set salt radius array
583  xs[ jj ] = sa_data.radius( jj );
584  }
585 DbgLv(2) << "SaltD: Nx" << Nx << "xs sme" << xs[0] << xs[1] << xs[2]
586  << xs[Nx/2-1] << xs[Nx/2] << xs[Nx-2+1] << xs[Nx-3] << xs[Nx-2] << xs[Nx-1];
587 };
588 
589 US_LammAstfvm::SaltData::~SaltData()
590 {
591  //delete [] xs;
592  //delete [] Cs0;
593  //delete [] Cs1;
594  xsVec .clear();
595  Cs0Vec.clear();
596  Cs1Vec.clear();
597 };
598 
599 void US_LammAstfvm::SaltData::initSalt()
600 {
601  t0 = sa_data.scanData[ 0 ].seconds; // times of 1st 2 salt scans
602  t1 = sa_data.scanData[ 1 ].seconds;
603  scn = 2; // index to next scan to use
604  Nt = sa_data.scanCount() - 2; // scan count less two used here
605 
606  for ( int j = 0; j < Nx; j++ )
607  { // get 1st two salt arrays from 1st two salt scans
608  Cs0[ j ] = sa_data.value( 0, j );
609  Cs1[ j ] = sa_data.value( 1, j );
610  }
611 int k=Nx/2;
612 int n=Nx-1;
613 DbgLv(2) << "initSalt: t0 t1" << t0 << t1;
614 DbgLv(2) << "initSalt: xs Cs0 Cs1 j" << xs[0] << Cs0[0] << Cs1[0] << 0;
615 DbgLv(2) << "initSalt: xs Cs0 Cs1 j" << xs[k] << Cs0[k] << Cs1[k] << k;
616 DbgLv(2) << "initSalt: xs Cs0 Cs1 j" << xs[n] << Cs0[n] << Cs1[n] << n;
617 }
618 
619 void US_LammAstfvm::SaltData::InterpolateCSalt( int N, double *x, double t,
620  double *Csalt )
621 {
622  double* tmp;
623 DbgLv(2) << "SaltD:ntrp: N t t1 Nt Nx" << N << t << t1 << Nt << Nx
624  << "Cs0N Cs1N" << Cs0[Nx-1] << Cs1[Nx-1];
625 
626  while ( ( t1 < t ) && ( Nt > 0 ) )
627  { // walk down salt scans until we are straddling desired time value
628  t0 = t1;
629  tmp = Cs0;
630  Cs0 = Cs1;
631  Cs1 = tmp; // swap Cs0 and Cs1
632 
633  t1 = sa_data.scanData[ scn ].seconds;
634 
635  for ( int j = 0; j < Nx; j++ )
636  Cs1[ j ] = sa_data.value( scn, j );
637 
638  Nt --; // Nt = time level left
639  scn++;
640 DbgLv(3) << "SaltD:ntrp: 0 t 1" << t0 << t << t1 << " N s" << Nt << scn;
641  }
642 DbgLv(2) << "SaltD:ntrp: t0 t t1" << t0 << t << t1 << " Nt scn" << Nt << scn;
643 
644  // interpolate between t0 and t1
645  double et1 = ( t - t0 ) / ( t1 - t0 );
646  et1 = ( et1 > 1.0 ) ? 1.0 : et1;
647  et1 = ( et1 < 0.0 ) ? 0.0 : et1;
648  double et0 = 1.0 - et1;
649 
650  // interpolate between xs[k-1] and xs[k]
651  int k = 1;
652  int Lx = Nx - 1;
653 int knan=0;
654 
655  for ( int j = 0; j < N; j++ ) // loop for all x[m]
656  {
657  double xj = x[ j ];
658  while ( xj > xs[ k ] && k < Lx ) k++; // radial point
659 
660  // linear interpolation
661  int m = k - 1;
662  double xik = ( xj - xs[ m ] ) / ( xs[ k ] - xs[ m ] );
663  xik = ( xik > 1.0 ) ? 1.0 : xik;
664  xik = ( xik < 0.0 ) ? 0.0 : xik;
665  double xim = 1.0 - xik;
666  // interpolate linearly in both time and radius
667  Csalt[ j ] = et0 * ( xim * Cs0[ m ] + xik * Cs0[ k ] )
668  + et1 * ( xim * Cs1[ m ] + xik * Cs1[ k ] );
669 //*DEBUG
670 double csj=Csalt[j];
671 if(csj!=csj) {
672 knan++;
673 if(knan<20)
674 DbgLv(2) << "SaltD:ntrp: j,k,xj,xsm,xsk" << j << k << xj << xs[m] << xs[k]
675  << "xim xik" << xim << xik << "Csaltj" << Csalt[j] << csj;
676 }
677 //*DEBUG
678  }
679 
680 DbgLv(2) << "SaltD:ntrp: scn" << scn << "N Nx" << N << Nx
681  << "Csalt0 CsaltM CsaltN" << Csalt[0] << Csalt[N/2] << Csalt[N-1];
682 DbgLv(2) << "SaltD:ntrp: Cs00 Cs0M Cs0N" << Cs0[0] << Cs0[Nx/2] << Cs0[Nx-1];
683 DbgLv(2) << "SaltD:ntrp: Cs10 Cs1M Cs1N" << Cs1[0] << Cs1[Nx/2] << Cs1[Nx-1];
684 }
685 
686 
687 // construct Lamm solver for finite volume method
688 US_LammAstfvm::US_LammAstfvm( US_Model& rmodel,
689  US_SimulationParameters& rsimparms,
690  QObject* parent /*=0*/ )
691  : QObject( parent ), model( rmodel ), simparams( rsimparms )
692 {
693  comp_x = 0; // initial model component index
694 
695  dbg_level = US_Settings::us_debug();
696  stopFlag = false;
697  movieFlag = false;
698  double speed = simparams.speed_step[ 0 ].rotorspeed;
699  double* coefs = simparams.rotorcoeffs;
700  double stretch = coefs[ 0 ] * speed + coefs[ 1 ] * sq( speed );
701  param_m = simparams.meniscus;
702  param_b = simparams.bottom;
703  param_w2 = sq( speed * M_PI / 30.0 );
704 DbgLv(2) << "LFvm: m b w2 strtch" << param_m << param_b << param_w2 << stretch;
705  param_m -= stretch;
706  param_b += stretch;
707 DbgLv(2) << "LFvm: m b" << param_m << param_b;
708 
709  MeshSpeedFactor = 1; // default mesh moving option
710  MeshRefineOpt = 1; // default mesh refinement option
711 
712  NonIdealCaseNo = 0; // default case: ideal with constant s, D
713 
714  sigma = 0.; // default: s=s_0
715  delta = 0.; // default: D=D_0
716  density = DENS_20W;
717  compressib = 0.0;
718  d_coeff[ 0 ] = DENS_20W;
719  v_coeff[ 0 ] = VISC_20W;
720  for ( int ii = 1; ii < 6; ii++ )
721  {
722  d_coeff[ ii ] = 0.0;
723  v_coeff[ ii ] = 0.0;
724  }
725 
726  //err_tol = 1.0e-4;
727  err_tol = 1.0e-5;
728 }
729 
730 // destroy
731 US_LammAstfvm::~US_LammAstfvm()
732 {
733  if ( NonIdealCaseNo == 2 ) delete saltdata;
734 
735  return;
736 }
737 
738 // primary method to calculate solutions for all species
739 int US_LammAstfvm::calculate( US_DataIO::RawData& sim_data )
740 {
741  auc_data = &sim_data;
742 
743  // use given data to create form for internal data; zero initial concs.
744  load_mfem_data( sim_data, af_data, true );
745 
746  // set up to report progress to any listener (e.g., us_astfem_sim)
747  int nsteps = af_data.scan.size() * model.components.size();
748 
749  if ( model.coSedSolute >= 0 )
750  { // for co-sedimenting, reduce total steps by one scan
751  nsteps -= af_data.scan.size();
752  }
753 
754 #ifndef NO_DB
755  emit calc_start( nsteps );
756  qApp->processEvents();
757 #endif
758 
759  // update concentrations for each model component
760  for ( int ii = 0; ii < model.components.size(); ii++ )
761  {
762  int rc = solve_component( ii );
763 
764  if ( rc != 0 )
765  return rc;
766  }
767 
768 #ifndef NO_DB
769  emit calc_done();
770  qApp->processEvents();
771 #endif
772 
773  // populate user's data set from calculated simulation
774  store_mfem_data( sim_data, af_data );
775 
776  return 0;
777 }
778 
779 // get a solution for a component and update concentrations
780 int US_LammAstfvm::solve_component( int compx )
781 {
782  comp_x = compx;
783 
784  param_s = model.components[ compx ].s;
785  param_D = model.components[ compx ].D;
786 
787  double t0 = 0.;
788  double t1 = 100.;
789  double* x0;
790  double* u0;
791  double* x1;
792  double* u1;
793  double* u1p0;
794  double* u1p;
795  double* dtmp;
796  double total_t = ( param_b - param_m ) * 2.0
797  / ( param_s * param_w2 * param_m );
798  double dt = log( param_b / param_m )
799  / ( param_w2 * param_s * 100.0 );
800 
801  int ntcc = (int)( total_t / dt ) + 1; // nbr. times in calculations
802  int jt = 0;
803  int nts = af_data.scan.size(); // nbr. output times (scans)
804  int kt = 0;
805  int ncs = af_data.radius.size(); // nbr. concentrations each scan
806  int N0;
807  int N1;
808  int N0u;
809  int N1u;
810  int istep = comp_x * nts;
811 
812  double solut_t = af_data.scan[ nts - 1 ].time; // true total time
813  int ntc = (int)( solut_t / dt ) + 1; // nbr. times in calculations
814 
815  QVector< double > conc0;
816  QVector< double > conc1;
817  QVector< double > rads;
818 
819 #ifndef NO_DB
820  emit comp_progress( compx + 1 );
821  qApp->processEvents();
822 #endif
823 
824  if ( nonIdealCaseNo() != 0 ) // set non-ideal case number
825  { // if multiple cases, abort
826  return 1;
827  }
828 
829  if ( NonIdealCaseNo == 3 )
830  { // compressibility: 8-fold smaller delta-t and greater time points
831  dt /= 8.0;
832  ntc = (int)( solut_t / dt ) + 1;
833  }
834 
835 DbgLv(2) << "LAsc: CX=" << comp_x
836  << " ntcc ntc nts ncs nicase" << ntcc << ntc << nts << ncs << NonIdealCaseNo;
837 DbgLv(2) << "LAsc: tot_t dt sol_t" << total_t << dt << solut_t;
838 DbgLv(2) << "LAsc: m b s w2" << param_m << param_b << param_s << param_w2;
839 
840  if ( NonIdealCaseNo == 2 )
841  { // co-sedimenting
842  if ( comp_x == model.coSedSolute )
843  { // if this component is the salt, skip solving for it
844  if ( comp_x == 0 )
845  {
846 DbgLv(2) << "NonIdeal2: new saltdata comp_x" << comp_x;
847  saltdata = new SaltData( model, simparams, auc_data );
848  vbar_salt = model.components[ 0 ].vbar20;
849  }
850 
851  //saltdata->initSalt();
852  return 0;
853  }
854 
855  else if ( compx > model.coSedSolute )
856  { // if beyond the salt component, adjust step count
857  istep -= nts;
858  }
859  }
860 
861  conc0.resize( ncs );
862  conc1.resize( ncs );
863  rads. resize( ncs );
864 
865  Mesh *msh = new Mesh( param_m, param_b, 100, 0 );
866 
867  msh->InitMesh( param_s, param_D, param_w2 );
868  int mropt = 1; // mesh refine option;
869 
870  // make settings based on non-ideal case type
871  if ( NonIdealCaseNo == 1 ) // concentration-dependent
872  {
873  SetNonIdealCase_1( model.components[ comp_x ].sigma,
874  model.components[ comp_x ].delta );
875 DbgLv(2) << "LAsc: sigma delta" << model.components[comp_x].sigma
876  << model.components[comp_x].delta << " comp_x" << comp_x;
877  }
878 
879  else if ( NonIdealCaseNo == 2 ) // co-sedimenting
880  {
881  SetNonIdealCase_2( );
882  }
883 
884  else if ( NonIdealCaseNo == 3 ) // compressibility
885  {
886  SetNonIdealCase_3( mropt, err_tol );
887  }
888 
889  else
890  {
891  NonIdealCaseNo = 0;
892  }
893 
894  SetMeshRefineOpt( mropt ); // mesh refine option
895  SetMeshSpeedFactor( 1.0 ); // mesh speed factor
896 
897 
898  // get initial concentration for this component
899  double sig_conc = model.components[ comp_x ].signal_concentration;
900 
901 QTime timer;
902 int ktime1=0;
903 int ktime2=0;
904 int ktime3=0;
905 int ktime4=0;
906 int ktime5=0;
907 int ktime6=0;
908 timer.start();
909  // initialization
910  N0 = msh->Nv;
911  N0u = N0 + N0 - 1;
912  x0 = new double [ N0 ];
913  u0 = new double [ N0u ];
914  N1 = N0;
915  N1u = N0u;
916  x1 = new double [ N1 ];
917  u1 = new double [ N1u ];
918 
919  for ( int jj = 0; jj < N0; jj++ )
920  { // initialize X and U values
921  int kk = jj + jj;
922 
923  x0[ jj ] = msh->x[ jj ]; // r value
924  x1[ jj ] = x0[ jj ];
925  u0[ kk ] = msh->x[ jj ] * sig_conc; // C*r value
926  u1[ kk ] = u0[ kk ];
927  }
928 
929  for ( int kk = 1; kk < N0u - 1; kk+=2 )
930  { // fill in mid-point U values
931  u0[ kk ] = ( u0[ kk - 1 ] + u0[ kk + 1 ] ) * 0.5;
932  u1[ kk ] = u0[ kk ];
933  }
934 DbgLv(2) << "LAsc: u0 0,1,2...,N" << u0[0] << u0[1] << u0[2]
935  << u0[N0u-3] << u0[N0u-2] << u0[N0u-1];
936 
937  for ( int jj = 0; jj < ncs; jj++ )
938  { // get output radius vector
939  rads[ jj ] = af_data.radius[ jj ];
940  }
941 DbgLv(2) << "LAsc: r0 rn ncs rsiz" << rads[0] << rads[ncs-1]
942  << ncs << af_data.radius.size() << rads.size();
943 
944  const double u0max = 1e+100;
945  const double u0min = -u0max;
946 
947  for ( int jj = 0; jj < N0u; jj++ )
948  {
949  u0[ jj ] = qMin( u0max, qMax( u0min, u0[ jj ] ) );
950  u1[ jj ] = qMin( u0max, qMax( u0min, u1[ jj ] ) );
951  }
952 DbgLv(2) << "LAsc: u0 0,1,2...,N" << u0[0] << u0[1] << u0[2]
953  << u0[N0u-3] << u0[N0u-2] << u0[N0u-1];
954 
955 #ifndef NO_DB
956  int ktinc = 5; // signal progress every 5th scan
957  if ( nts > 10000 )
958  ktinc = qMax( 2, ( ( nts + 5000 ) / 10000 ) );
959  if ( nts < 100 )
960  ktinc = 1;
961 DbgLv(2) << "LAsc: nts ktinc" << nts << ktinc;
962 #endif
963  double ts;
964  double u_ttl;
965 
966 
967  QDir dir;
968  QString tmpDir = US_Settings::tmpDir();
969  if ( ! dir.exists( tmpDir ) ) dir.mkpath( tmpDir );
970 
971  QFile ftto( tmpDir + "/tt0-ufvm" );
972  if ( dbg_level > 0 )
973  {
974  if ( ! ftto.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
975  qDebug() << "*ERROR* Unable to open tt0-ufvm";
976  }
977  QTextStream tso( &ftto );
978 //ntc=ntcc;
979  if ( dbg_level > 0 )
980  tso << ntc << "\n";
981 
982  // main loop for time
983  for ( jt = 0, kt = 0; jt < ntc; jt++ )
984  {
985 timer.restart();
986  t0 = dt * (double)jt;
987  t1 = t0 + dt;
988  ts = af_data.scan[ kt ].time; // time at output scan
989  N0u = N0 + N0 - 1;
990  if ( dbg_level > 0 && ( ( jt / 10 ) * 10 ) == jt )
991  {
992  u_ttl = IntQs( x0, u0, 0, -1, N0-2, 1 );
993 DbgLv(2) << "LAsc: jt,kt,t0=" << jt << kt << t0 << " Nv=" << N0
994  << "u_ttl=" << u_ttl;
995 DbgLv(2) << "LAsc: u0 0,1,2...,N" << u0[0] << u0[1] << u0[2]
996  << u0[N0u-3] << u0[N0u-2] << u0[N0u-1];
997  tso << QString().sprintf( "%12.5e %d %12.5e\n", t0, N0, u_ttl );
998  for ( int j=0; j<N0; j++ )
999  tso << QString().sprintf( "%10.6e \n", x0[j] );
1000  tso << QString().sprintf( "\n" );
1001  for ( int j=0; j<N0u; j++ )
1002  tso << QString().sprintf( "%15.7e \n", u0[j] );
1003  tso << QString().sprintf( "\n\n" );
1004  }
1005 ktime1+=timer.restart();
1006 
1007  u1p0 = new double [ N0u ];
1008  LammStepSedDiff_P( t0, dt, N0-1, x0, u0, u1p0 );
1009 ktime2+=timer.restart();
1010 
1011  if ( MeshRefineOpt == 1 )
1012  {
1013 
1014  msh->RefineMesh( u0, u1p0, err_tol );
1015 ktime3+=timer.restart();
1016 
1017  N1 = msh->Nv;
1018  N1u = N1 + N1 - 1;
1019  u1p = new double [ N1u ];
1020 
1021  delete [] x1;
1022  x1 = new double [ N1 ];
1023 
1024  for ( int jj = 0; jj < N1; jj++ )
1025  x1[ jj ] = msh->x[ jj ];
1026 
1027  ProjectQ( N0-1, x0, u1p0, N1-1, x1, u1p );
1028 
1029  delete [] u1;
1030  u1 = new double [ N1u ];
1031 
1032 ktime4+=timer.restart();
1033  LammStepSedDiff_C( t0, dt, N0-1, x0, u0, N1-1, x1, u1p, u1 );
1034 
1035  delete [] u1p;
1036  }
1037 
1038  else
1039  {
1040 ktime4+=timer.restart();
1041  LammStepSedDiff_C( t0, dt, N0-1, x0, u0, N1-1, x1, u1p0, u1 );
1042  }
1043 
1044 ktime5+=timer.restart();
1045  // see if current scan is between calculated times; output scan if so
1046 
1047  if ( ts >= t0 && ts <= t1 )
1048  { // interpolate concentrations quadratically in x; linearly in time
1049  double f0 = ( t1 - ts ) / ( t1 - t0 ); // fraction of conc0
1050  double f1 = ( ts - t0 ) / ( t1 - t0 ); // fraction of conc1
1051 
1052 DbgLv(2) << "LAsc: call qI t0 ts t1" << t0 << ts << t1;
1053  // do quadratic interpolation to fill out concentrations at time t0
1054  quadInterpolate( x0, u0, N0, rads, conc0 );
1055 
1056  // do quadratic interpolation to fill out concentrations at time t1
1057  quadInterpolate( x1, u1, N1, rads, conc1 );
1058 DbgLv(2) << "LAsc: x0[0] x0[H] x0[N]" << x0[0] << x0[N0/2] << x0[N0-1];
1059 DbgLv(2) << "LAsc: x1[0] x1[H] x1[N]" << x1[0] << x1[N1/2] << x1[N1-1];
1060 DbgLv(2) << "LAsc: r[0] r[H] r[N]" << rads[0] << rads[ncs/2] << rads[ncs-1];
1061 DbgLv(2) << "LAsc: u0[0] u0[H] u0[N]" << u0[0] << u0[N0u/2] << u0[N1u-1];
1062 DbgLv(2) << "LAsc: u1[0] u1[H] u1[N]" << u1[0] << u1[N1u/2] << u1[N1u-1];
1063 DbgLv(2) << "LAsc: c0[0] c0[H] c0[N]"
1064  << conc0[0] << conc0[ncs/2] << conc0[ncs-1];
1065 DbgLv(2) << "LAsc: c1[0] c1[H] c1[N]"
1066  << conc1[0] << conc1[ncs/2] << conc1[ncs-1];
1067  double utt0 = IntQs( x0, u0, 0, -1, N0-2, 1 );
1068  double utt1 = IntQs( x1, u1, 0, -1, N1-2, 1 );
1069 DbgLv(2) << "LAsc: utt0 utt1" << utt0 << utt1;
1070 DbgLv(2) << "LAsc:stopFlag" << stopFlag;
1071 
1072  double cmax = 0.0;
1073  double rmax = 0.0;
1074 
1075 //f0=1.0; f1=0.0;
1076  for ( int jj = 0; jj < ncs; jj++ )
1077  { // update concentration vector with linear interpolation for time
1078  af_data.scan[ kt ].conc[ jj ] += ( conc0[ jj ] * f0 +
1079  conc1[ jj ] * f1 );
1080  double Cm = af_data.scan[ kt ].conc[ jj ];
1081  if ( Cm > cmax )
1082  {
1083  cmax = Cm;
1084  rmax = af_data.radius[ jj ];
1085  }
1086  }
1087 DbgLv(2) << "LAsc: t=" << ts << "Cmax=" << cmax << " r=" << rmax;
1088 DbgLv(2) << "LAsc: co[0] co[H] co[N] kt" << af_data.scan[kt].conc[0]
1089  << af_data.scan[kt].conc[ncs/2] << af_data.scan[kt].conc[ncs-1] << kt;
1090 
1091  istep++; // bump progress step
1092 
1093 #ifndef NO_DB
1094  if ( ( ( kt / ktinc ) * ktinc ) == kt || ( kt + 1 ) == nts )
1095  { // signal progress at every "ktinc'th" scan or final one
1096  emit calc_progress( istep );
1097 DbgLv(2) << "LAsc: istep" << istep;
1098  qApp->processEvents();
1099  }
1100 #endif
1101 #ifndef NO_DB
1102 // if ( movieFlag &&
1103 // ( ( ( kt / ktinc ) * ktinc ) == kt || ( kt + 1 ) == nts ) )
1104  if ( movieFlag )
1105  {
1106  emit new_scan( &af_data.radius, af_data.scan[ kt ].conc.data() );
1107  emit new_time( t0 );
1108  qApp->processEvents();
1109  }
1110 #endif
1111 
1112  kt++; // bump output time(scan) index
1113 
1114  qApp->processEvents();
1115  if ( stopFlag ) break;
1116  }
1117 
1118  delete [] u1p0;
1119 
1120  if ( kt >= nts )
1121  break; // if all scans updated, we are done
1122 
1123  qApp->processEvents();
1124  if ( stopFlag ) break;
1125  // switch x,u arrays for next iteration
1126  N0 = N1;
1127  N0u = N1u;
1128  dtmp = x0;
1129  x0 = x1;
1130  x1 = dtmp;
1131  dtmp = u0;
1132  u0 = u1;
1133  u1 = dtmp;
1134  }
1135 
1136  if ( dbg_level > 0 )
1137  ftto.close();
1138 
1139  if ( dbg_level > 0 )
1140  { // calculate and print the integral of scan curves
1141  double cimn = 9e+14;
1142  double cimx = 0.0;
1143  double ciav = 0.0;
1144  double dltr = ( af_data.radius[ 1 ] - af_data.radius[ 0 ] ) * 0.5;
1145 
1146  for ( int ii = 0; ii < af_data.scan.size(); ii++ )
1147  {
1148  double csum = 0.0;
1149  double pval = af_data.scan[ ii ].conc[ 0 ];
1150 
1151  for ( int jj = 1; jj < af_data.scan[ ii ].conc.size(); jj++ )
1152  {
1153  double cval = af_data.scan[ ii ].conc[ jj ];
1154  csum += ( ( cval + pval ) * dltr );
1155  pval = cval;
1156 //if ( ii < 19 && ( (jj/100)*100 == jj || (jj+5)>nconc ) )
1157 // DbgLv(3) << " jj cval dltr csum" << jj << cval << dltr << csum;
1158  }
1159 
1160  DbgLv(2) << "Scan" << ii + 1 << " Integral" << csum;
1161  cimn = ( cimn < csum ) ? cimn : csum;
1162  cimx = ( cimx > csum ) ? cimx : csum;
1163  ciav += csum;
1164  }
1165 
1166  ciav /= (double)af_data.scan.size();
1167  double cidf = cimx - cimn;
1168  double cidp = (double)( qRound( 10000.0 * cidf / ciav ) ) / 100.0;
1169  DbgLv(2) << " Integral Min Max Mean" << cimn << cimx << ciav;
1170  DbgLv(2) << " ( range of" << cidf << "=" << cidp << " percent of mean )";
1171  }
1172 
1173  delete [] x0; // clean up
1174  delete [] u0;
1175  delete [] x1;
1176  delete [] u1;
1177  delete msh;
1178 ktime6+=timer.elapsed();
1179 DbgLv(2) << "compx" << comp_x << "times 1-6"
1180  << ktime1 << ktime2 << ktime3 << ktime4 << ktime5 << ktime6;
1181  return 0;
1182 }
1183 
1184 void US_LammAstfvm::set_buffer( US_Buffer buffer )
1185 {
1186  density = buffer.density; // for compressibility
1187  compressib = buffer.compressibility;
1188 
1189  buffer.compositeCoeffs( d_coeff, v_coeff );
1190 DbgLv(2) << "buff d_coeff" << d_coeff[0] << d_coeff[1] << d_coeff[2]
1191  << d_coeff[3] << d_coeff[4] << d_coeff[5];
1192 DbgLv(2) << "buff v_coeff" << v_coeff[0] << v_coeff[1] << v_coeff[2]
1193  << v_coeff[3] << v_coeff[4] << v_coeff[5];
1194 }
1195 
1196 void US_LammAstfvm::SetNonIdealCase_1( double sigma_k, double delta_k )
1197 {
1198  sigma = sigma_k; // for concentration dependency
1199  delta = delta_k;
1200 }
1201 
1202 void US_LammAstfvm::SetNonIdealCase_2( )
1203 {
1204  if ( comp_x == 0 )
1205  {
1206 DbgLv(2) << "NonIdeal2: create saltdata";
1207  saltdata = new SaltData( model, simparams, auc_data );
1208  vbar_salt = model.components[ 0 ].vbar20;
1209  }
1210 
1211 DbgLv(2) << "NonIdeal2: initSalt comp_x" << comp_x;
1212  saltdata->initSalt();
1213 }
1214 
1215 void US_LammAstfvm::SetNonIdealCase_3( int& mropt, double& err_tol )
1216 {
1217  mropt = 0;
1218  err_tol = 1.0e-5;
1219 }
1220 
1221 void US_LammAstfvm::SetMeshSpeedFactor( double speed )
1222 {
1223  MeshSpeedFactor = speed;
1224 }
1225 
1226 void US_LammAstfvm::SetMeshRefineOpt( int Opt )
1227 {
1228  MeshRefineOpt = Opt;
1229 }
1230 
1232 //
1233 // LammStepSedDiff_P:
1234 // prediction step of solving Lamm equation (sedimentation-diffusion only)
1235 //
1236 // Given solution (x0, u0) at t, to find solution (x1, u1) at time t+dt
1237 // use u1p = u0
1238 //
1239 // M0 = number of elems in x0
1240 // u0 = piecewise quadratic solution at t on mesh x0
1241 // u1 = piecewise quadratic solution at t+dt on mesh x1
1242 //
1244 void US_LammAstfvm::LammStepSedDiff_P( double t, double dt, int M0,
1245  double *x0, double *u0, double *u1 )
1246 {
1247  LammStepSedDiff_C( t, dt, M0, x0, u0, M0, x0, u0, u1 );
1248 }
1249 
1250 
1252 //
1253 // LammStepSedDiff_C:
1254 // Correction step of solving Lamm equation (sedimentation-diffusion only)
1255 //
1256 // Given solution (x0, u0) at t, and estimate (x0, u_est) of
1257 // solution at t+dt, and mesh x1, to find solution (x1, u1)
1258 // at time t+dt
1259 //
1260 // M0, M1 = number of elems in x0, x1
1261 // u0 = piecewise quadratic solution at t on mesh x0
1262 // u1 = piecewise quadratic solution at t+dt on mesh x1
1263 // u1p = piecewise quadratic estimated solution at t+dt on mesh x1
1264 //
1266 void US_LammAstfvm::LammStepSedDiff_C( double t, double dt, int M0,
1267  double *x0, double *u0, int M1, double *x1, double *u1p, double *u1 )
1268 {
1269  int Ng = 2 * M1; // number of x_star points
1270  int* ke = new int [ Ng ];
1271  double* MemDouble = new double [ 12 * Ng + 15 ];
1272  double* flux_p[ 3 ];
1273 
1274  double dt2 = dt * 0.5;
1275  double* xt = MemDouble;
1276  double* xi = xt + Ng;
1277  double* xg0 = xi + Ng;
1278  double* xg1 = xg0 + Ng;
1279  double* ug0 = xg1 + Ng;
1280  double* ug1 = ug0 + Ng;
1281  double* Sv = ug1 + Ng;
1282  double* Dv = Sv + Ng;
1283  double* flux_u = Dv + Ng;
1284  flux_p[ 0 ] = flux_u + Ng;
1285  flux_p[ 1 ] = flux_p[ 0 ] + Ng;
1286  flux_p[ 2 ] = flux_p[ 1 ] + Ng;
1287  double* phi = flux_p[ 2 ] + Ng;
1288  double* phiL = phi + 3;
1289  double* phiR = phiL + 6;
1290 QTime timer;
1291 static int ktim1=0;
1292 static int ktim2=0;
1293 static int ktim3=0;
1294 static int ktim4=0;
1295 static int ktim5=0;
1296 static int ktim6=0;
1297 static int ktim7=0;
1298 static int ktim8=0;
1299 timer.start();
1300 
1301  // calculate Sv, Dv at t+dt on xg=(xl, xr)
1302  for ( int j = 0; j < Ng; j += 2 )
1303  {
1304  int j2 = j / 2;
1305  xg1[ j ] = x1[ j2 ] * 0.75 + x1[ j2 + 1 ] * 0.25; // xl
1306  xg1[ j + 1 ] = x1[ j2 ] * 0.25 + x1[ j2 + 1 ] * 0.75; // xr
1307  ug1[ j ] = ( 3. * u1p[ j ] + 6. * u1p[ j + 1 ] - u1p[ j + 2 ] ) / 8.;
1308  ug1[ j + 1 ] = ( 3. * u1p[ j + 2 ] + 6. * u1p[ j + 1 ] - u1p[ j ] ) / 8.;
1309  }
1310 
1311  AdjustSD( t + dt, Ng, xg1, ug1, Sv, Dv );
1312 ktim1+=timer.restart();
1313 DbgLv(2) << " xg1 0 1 M Nm N" << xg1[0] << xg1[1] << xg1[Ng/2]
1314  << xg1[Ng-2] << xg1[Ng-1];
1315 DbgLv(2) << " Sv 0 1 M Nm N" << Sv[0] << Sv[1] << Sv[Ng/2]
1316  << Sv[Ng-2] << Sv[Ng-1];
1317 
1318  // determine xg0=(xls, xrs)
1319  for ( int j = 0; j < Ng; j++ )
1320  {
1321  double sw2 = Sv[ j ] * param_w2;
1322  xg0[ j ] = xg1[ j ] - dt * MeshSpeedFactor * sw2 * xg1[ j ]
1323  * exp( -qAbs( sw2 ) * dt / 2 );
1324 
1325  xg0[ j ] = qMax( param_m, qMin( param_b, xg0[ j ] ) );
1326 
1327  xt[ j ] = ( xg1[ j ] - xg0[ j ] ) / dt;
1328  }
1329 DbgLv(2) << " xg0 0 1 M Nm N" << xg0[0] << xg0[1] << xg0[Ng/2]
1330  << xg0[Ng-2] << xg0[Ng-1] << "Ng" << Ng;
1331 
1332  // redistribute xgs so that in between [m,b] and in increasing order
1333  double bl = param_m;
1334 
1335  for ( int j = 0; j < Ng; j++ )
1336  {
1337  int cnt = 1;
1338 
1339  while ( xg0[ j ] < bl && ( j + 1 ) < Ng )
1340  {
1341  j++;
1342  cnt++;
1343  }
1344 
1345  double br = qMin( xg0[ j ], param_b );
1346 
1347  for ( int jm = 0; jm < cnt; jm++ )
1348  {
1349  xg0[ j - jm ] = br - (double)jm / (double)cnt * ( br - bl );
1350  xt[ j - jm ] = ( xg1[ j - jm ] - xg0[ j - jm ] ) / dt;
1351  }
1352 
1353  bl = br;
1354  }
1355 DbgLv(2) << " xg0 0 1 M Nm N" << xg0[0] << xg0[1] << xg0[Ng/2]
1356  << xg0[Ng-2] << xg0[Ng-1] << "Ng" << Ng;
1357 ktim2+=timer.restart();
1358 
1359  // calculate Flux(phi, t+dt) at all xg1
1360  fun_dphi( -0.5, phiL ); // basis at xi=-1/2
1361  fun_dphi( 0.5, phiR );
1362 
1363  for ( int j = 0; j < Ng; j++ )
1364  {
1365  int j2 = j / 2;
1366  double h = 0.5 * ( x1[ j2 + 1 ] - x1[ j2 ] );
1367 
1368  for ( int jm = 0; jm < 3; jm++ ) // at xl
1369  {
1370  flux_p[ jm ][ j ] = ( xt[ j ] - Sv[ j ] * param_w2 * xg1[ j ]
1371  - Dv[ j ] / xg1[ j ] ) * phiL[ jm ]
1372  + Dv[ j ] * phiL[ jm + 3 ] / h;
1373  }
1374  j++;
1375 
1376  for ( int jm = 0; jm < 3; jm++ ) // at xr
1377  {
1378  flux_p[ jm ][ j ] = ( xt[ j ] - Sv[ j ] * param_w2 * xg1[ j ]
1379  - Dv[ j ] / xg1[ j ] ) * phiR[ jm ]
1380  + Dv[ j ] * phiR[ jm + 3 ] / h;
1381  }
1382  }
1383 
1384 
1385  // calculate Sv, Dv at (xg0, t)
1386 
1387  LocateStar( M0 + 1, x0, Ng, xg0, ke, xi ); // position of xg0 on mesh x0
1388 
1389  for ( int j = 0; j < Ng; j++ )
1390  {
1391  fun_phi( xi[ j ], phi );
1392 
1393  int j2 = 2 * ke[ j ];
1394  ug0[ j ] = u0[ j2 ] * phi[ 0 ]
1395  + u0[ j2 + 1 ] * phi[ 1 ]
1396  + u0[ j2 + 2 ] * phi[ 2 ];
1397  }
1398 
1399  // calculate s, D at xg0 on time t
1400 ktim3+=timer.restart();
1401  AdjustSD( t, Ng, xg0, ug0, Sv, Dv );
1402 ktim4+=timer.restart();
1403 
1404  // calculate Flux(u0,t) at all xg0
1405  // (i) Compute ux at nodes as average of Du from left and right
1406 
1407  double* ux = new double [ M0 + 1 ]; // D_x(u0) at all x0
1408 
1409  for ( int j = 1; j < M0; j++ ) // internal nodes
1410  {
1411  int j2 = 2 * j;
1412  ux[ j ] = ( ( u0[ j2 - 2 ] - 4. * u0[ j2 - 1 ] + 3. * u0[ j2 ] )
1413  / ( x0[ j ] - x0[ j - 1 ] )
1414  - ( 3. * u0[ j2 ] - 4. * u0[ j2 + 1 ] + u0[ j2 + 2 ] )
1415  / ( x0[ j + 1 ] - x0[ j ] ) ) / 2.;
1416  }
1417 
1418  int j2 = 2 * M0;
1419  ux[ 0 ] = -( 3. * u0[ 0 ] - 4. * u0[ 1 ] + u0[ 2 ] )
1420  / ( x0[ 1 ] - x0[ 0 ] ) / 2.;
1421  ux[ M0 ] = ( u0[ j2 - 2 ] - 4. * u0[ j2 - 1 ] + 3. * u0[ j2 ] )
1422  / ( x0[ M0 ] - x0[ M0 - 1 ] ) / 2.;
1423 
1424  // (ii) flux(u0,t) at all xg0
1425  for ( int j = 0; j < Ng; j++ )
1426  {
1427  double wt = ( 1. - xi[ j ] ) / 2.;
1428  double uxs = ux[ ke[ j ] ] * wt
1429  + ux[ ke[ j ] + 1 ] * ( 1. - wt ); // Du0 at xg0
1430 
1431  if ( ( xg0[ j ] <= ( param_m + 1.e-14 ) ) ||
1432  ( xg0[ j ] >= ( param_b - 1.e-14 ) ) )
1433  {
1434  flux_u[ j ] = 0.;
1435  }
1436 
1437  else
1438  {
1439  flux_u[ j ] = -( Sv[ j ] * param_w2 * xg0[ j ] + Dv[ j ] / xg0[ j ] )
1440  * ug0[ j ] + Dv[ j ] * uxs + xt[ j ] * ug0[ j ];
1441  }
1442  }
1443 
1444  delete [] ux;
1445 
1446  //
1447  // assemble the linear system of equations
1448  //
1449  double** Mtx = new double* [ Ng + 1 ];
1450  double* rhs = new double [ Ng + 1 ];
1451  for ( int i = 0; i <= Ng; i++ )
1452  Mtx[ i ] = new double [ 5 ];
1453 
1454 ktim5+=timer.restart();
1455  // Assemble the coefficient matrix
1456  for ( int i = 1; i < Ng; i += 2 )
1457  {
1458  int k = ( i - 1 ) / 2;
1459  double h = 0.5 * ( x1[ k + 1 ] - x1[ k ] );
1460  Mtx[ i ][ 0 ] = 0.;
1461  Mtx[ i ][ 1 ] = h / 24.
1462  + dt2 * ( flux_p[ 0 ][ i - 1 ] - flux_p[ 0 ][ i ] );
1463  Mtx[ i ][ 2 ] = h * 22. / 24.
1464  + dt2 * ( flux_p[ 1 ][ i - 1 ] - flux_p[ 1 ][ i ] );
1465  Mtx[ i ][ 3 ] = h / 24.
1466  + dt2 * ( flux_p[ 2 ][ i - 1 ] - flux_p[ 2 ][ i ] );
1467  Mtx[ i ][ 4 ] = 0.;
1468  }
1469 
1470  for ( int i = 2; i < Ng; i += 2 )
1471  {
1472  int k = i / 2;
1473 
1474  double h = 0.5 * ( x1[ k ] - x1[ k - 1 ] );
1475  Mtx[ i ][ 0 ] =-h / 24. + dt2 * flux_p[ 0 ][ i - 1 ];
1476  Mtx[ i ][ 1 ] = h * 5. / 24. + dt2 * flux_p[ 1 ][ i - 1 ];
1477  Mtx[ i ][ 2 ] = h * 8. / 24. + dt2 * flux_p[ 2 ][ i - 1 ];
1478 
1479  h = 0.5 * ( x1[ k + 1 ] - x1[ k ] );
1480  Mtx[ i ][ 2 ] += h * 8. / 24. - dt2 * flux_p[ 0 ][ i ];
1481  Mtx[ i ][ 3 ] = h * 5. / 24. - dt2 * flux_p[ 1 ][ i ];
1482  Mtx[ i ][ 4 ] =-h / 24. - dt2 * flux_p[ 2 ][ i ];
1483  }
1484 
1485  int i = 0;
1486  double h = 0.5 * ( x1[ 1 ] - x1[ 0 ] );
1487  Mtx[ i ][ 2 ] = h * 8. / 24. - dt2 * flux_p[ 0 ][ 0 ];
1488  Mtx[ i ][ 3 ] = h * 5. / 24. - dt2 * flux_p[ 1 ][ 0 ];
1489  Mtx[ i ][ 4 ] =-h / 24. - dt2 * flux_p[ 2 ][ 0 ];
1490 
1491  i = Ng;
1492  h = 0.5 * ( x1[ M1 ] - x1[ M1 - 1 ] );
1493  Mtx[ i ][ 0 ] =-h / 24. + dt2 * flux_p[ 0 ][ i - 1 ];
1494  Mtx[ i ][ 1 ] = h * 5. / 24. + dt2 * flux_p[ 1 ][ i - 1 ];
1495  Mtx[ i ][ 2 ] = h * 8. / 24. + dt2 * flux_p[ 2 ][ i - 1 ];
1496 
1497 ktim6+=timer.restart();
1498 
1499  // assemble the right hand side
1500  i = 0;
1501  rhs[ i ] = IntQs( x0, u0, 0, -1., ke[ i ], xi[ i ] )
1502  + dt2 * flux_u[ i ];
1503 
1504  for ( int i = 1; i < Ng; i++ )
1505  {
1506  rhs[i] = IntQs( x0, u0, ke[ i - 1 ], xi[ i - 1 ], ke[ i ], xi[ i ] )
1507  + dt2 * ( flux_u[ i ] - flux_u[i - 1 ] );
1508  }
1509 
1510  i = Ng;
1511  rhs[ i ] = IntQs( x0, u0, ke[ i - 1 ], xi[ i - 1 ], M0 - 1, 1. )
1512  + dt2 * ( - flux_u[ i - 1 ] );
1513 
1514 ktim7+=timer.restart();
1515  LsSolver53( Ng, Mtx, rhs, u1 );
1516 ktim8+=timer.restart();
1517 
1518  for ( int i = 0; i <= Ng; i++ )
1519  delete [] Mtx[ i ];
1520 
1521  delete [] Mtx;
1522  delete [] rhs;
1523 
1524 
1525  delete [] ke;
1526  delete [] MemDouble;
1527 DbgLv(2) << " Diff_C times 1-8" << ktim1 << ktim2 << ktim3 << ktim4
1528  << ktim5 << ktim6 << ktim7 << ktim8;
1529 }
1530 
1531 
1532 
1534 //
1535 // find in the mesh x0 the index of elem containing each xs,
1536 // and the xi-coordinates of xs in the elem
1537 //
1538 // Note: xs has to be in increasing order
1539 // N0 = number of points in x0
1540 // Ns = number of points in xs
1541 //
1543 
1544 void US_LammAstfvm::LocateStar( int N0, double *x0, int Ns, double *xs,
1545  int *ke, double *xi )
1546 {
1547  int eix = 1;
1548 
1549  for ( int j = 0; j < Ns; j++ )
1550  {
1551  while ( xs[ j ] > x0[ eix ] && eix < ( N0 - 1 ) ) eix ++;
1552 
1553  ke[ j ] = eix - 1;
1554  xi[ j ] = ( xs[ j ] - x0[ eix - 1 ] )
1555  / ( x0[ eix ] - x0[ eix - 1 ] ) * 2. - 1.;
1556  }
1557 }
1558 
1560 //
1561 // Find the adjusted s and D values according to time t, location x,
1562 // and concentration C = u/x
1563 // (note: the input is u=r*C)
1564 //
1566 void US_LammAstfvm::AdjustSD( double t, int Nv, double *x, double *u,
1567  double *s_adj, double *D_adj )
1568 {
1569  const double Tempt = 293.15; // temperature in K
1570  const double vbar_w = 0.72;
1571  const double rho_w = 0.998234; // density of water
1572  int jj;
1573  QVector< double > CsaltVec( Nv );
1574  double* Csalt;
1575  double Cm = 0.0;
1576  double rho;
1577  double visc;
1578  //double vbar = 0.251;
1579  //double vbar = 0.72; // 0.251;
1580  //double vbar = model.components[ 0 ].vbar20;
1581  double vbar = model.components[ comp_x ].vbar20;
1582 QTime timer;
1583 static int kst1=0;
1584 static int kst2=0;
1585 
1586  switch ( NonIdealCaseNo )
1587  {
1588  case 0: // ideal, s=s_0, D=D_0
1589 
1590  for ( jj = 0; jj < Nv; jj++ )
1591  {
1592  s_adj[ jj ] = param_s ;
1593  D_adj[ jj ] = param_D ;
1594  }
1595  break;
1596 
1597  case 1: // concentration dependent
1598  for ( jj = 0; jj < Nv; jj++ )
1599  {
1600  s_adj[ jj ] = param_s / ( 1. + sigma * u[ jj ] / x[ jj ] );
1601  D_adj[ jj ] = param_D / ( 1. + delta * u[ jj ] / x[ jj ] );
1602  }
1603  break;
1604 
1605  case 2: // co-sedimenting
1606  //** salt-protein
1607 timer.start();
1608 
1609  CsaltVec.resize( Nv );
1610  Csalt = CsaltVec.data();
1611 
1612 DbgLv(2) << "NonIdeal2: ntrp Salt";
1613  saltdata->InterpolateCSalt( Nv, x, t, Csalt ); // Csalt at (x, t)
1614 kst1+=timer.restart();
1615  {
1616 double rho0=0.0;
1617 double rhom=0.0;
1618 double rhoe=0.0;
1619 double vis0=0.0;
1620 double vism=0.0;
1621 double vise=0.0;
1622 double cms0=0.0;
1623 double cmsm=0.0;
1624 double cmse=0.0;
1625  double sA = param_s * 1.00194 / ( 1.0 - vbar * rho_w );
1626  double dA = param_D * Tempt * 1.00194 / 293.15;
1627  double Cmrt = 0.0;
1628  double Cmsq = 0.0;
1629  double Cmcu = 0.0;
1630  double Cmqu = 0.0;
1631 
1632  for ( jj = 0; jj < Nv; jj++ )
1633  {
1634  // The calculations below are a more efficient version of the
1635  // following original computations:
1636  //
1637  //rho = 0.998234 + Cm*( 12.68641e-2 + Cm*( 1.27445e-3 +
1638  // Cm*( -11.954e-4 + Cm*258.866e-6 ) ) ) + 6.e-6;
1639  //visc = 1.00194 - 19.4104e-3*sqrt(Cm) + Cm*( -4.07863e-2 +
1640  // Cm*( 11.5489e-3 + Cm*(-21.774e-4) ) ) - 0.00078;
1641  //s_adj[j] =(1-vbar*rho)*1.00194/((1-vbar_w*rho_w)*visc)*param_s;
1642  //D_adj[j] =(Tempt*1.00194)/(293.15*visc) * param_D;
1643 
1644  Cm = Csalt[ jj ]; // Salt concentration
1645  Cmrt = sqrt( Cm ); // and powers of it
1646  Cmsq = Cm * Cm;
1647  Cmcu = Cm * Cmsq;
1648  Cmqu = Cm * Cmcu;
1649 
1650  rho = d_coeff[ 0 ]
1651  + Cmrt * d_coeff[ 1 ]
1652  + Cm * d_coeff[ 2 ]
1653  + Cmsq * d_coeff[ 3 ]
1654  + Cmcu * d_coeff[ 4 ]
1655  + Cmqu * d_coeff[ 5 ];
1656 
1657  visc = v_coeff[ 0 ]
1658  + Cmrt * v_coeff[ 1 ]
1659  + Cm * v_coeff[ 2 ]
1660  + Cmsq * v_coeff[ 3 ]
1661  + Cmcu * v_coeff[ 4 ]
1662  + Cmqu * v_coeff[ 5 ];
1663  visc = qAbs( visc );
1664 
1665  s_adj[ jj ] = sA * ( 1.0 - vbar * rho ) / visc;
1666 
1667  D_adj[ jj ] = dA / visc;
1668  //D_adj[ jj ] = qAbs( dA / visc );
1669 
1670 if(jj==0){rho0=rho;vis0=visc;cms0=Cm;}
1671 if(jj==Nv/2){rhom=rho;vism=visc;cmsm=Cm;}
1672 if(jj==Nv-1){rhoe=rho;vise=visc;cmse=Cm;}
1673  }
1674 DbgLv(2) << "AdjSD: Csalt0 CsaltN Cm" << Csalt[0] << Csalt[Nv-1] << Cm;
1675 DbgLv(2) << "AdjSD: sadj 0 m n" << s_adj[0] << s_adj[Nv/2] << s_adj[Nv-1];
1676 DbgLv(2) << "AdjSD: Dadj 0 m n" << D_adj[0] << D_adj[Nv/2] << D_adj[Nv-1];
1677 DbgLv(2) << "AdjSD: rho 0,m,e" << rho0 << rhom << rhoe;
1678 DbgLv(2) << "AdjSD: visc 0,m,e" << vis0 << vism << vise;
1679 DbgLv(2) << "AdjSD: Cm 0,m,e" << cms0 << cmsm << cmse;
1680 DbgLv(2) << "AdjSD: vbar vbar_w rho_w" << vbar << vbar_w << rho_w;
1681 DbgLv(2) << "AdjSD: cmrt cm^2 cm^3 cm^4" << Cmrt << Cmsq << Cmcu << Cmqu;
1682 DbgLv(2) << "AdjSD: d_coeff[0] rho" << d_coeff[0] << rhoe;
1683 DbgLv(2) << "AdjSD: v_coeff[0] vis" << v_coeff[0] << vise;
1684  }
1685 
1686 kst2+=timer.restart();
1687 DbgLv(3) << "AdjSD: times 1 2" << kst1 << kst2;
1688  break;
1689 
1690  case 3: // compressibility
1691  {
1692  double phip = vbar; // apparent specific volume
1693  double alpha = 1.0;
1694  double factn = 0.5 * density * param_w2 * compressib;
1695  double msq = param_m * param_m;
1696  double sA = 1.0 - vbar * rho_w;
1697 
1698  for ( int jj = 0; jj < Nv; jj++ )
1699  {
1700  rho = density
1701  / ( 1.0 - factn * ( x[ jj ] * x[ jj ] - msq ) );
1702  double beta = ( 1.0 - phip * rho ) / sA;
1703  s_adj[ jj ] = param_s * alpha * beta;
1704  D_adj[ jj ] = param_D * alpha;
1705  }
1706 DbgLv(3) << "AdjSD: compr dens" << compressib << density;
1707 DbgLv(3) << "AdjSD: factn msq sa" << factn << msq << sA;
1708 DbgLv(3) << "AdjSD: sadj 0 m n" << s_adj[0] << s_adj[Nv/2] << s_adj[Nv-1];
1709 DbgLv(3) << "AdjSD: Dadj 0 m n" << D_adj[0] << D_adj[Nv/2] << D_adj[Nv-1];
1710  }
1711  break;
1712 
1713  default:
1714  qDebug( "invalid case number for non-ideal sedimentation" );
1715  break;
1716 
1717  } // switch
1718 
1719  return;
1720 }
1721 
1722 
1724 //
1725 // function [u1] = ProjectQ(x0, u0, x1)
1726 // Project piecewise quadratic solution u0 at x0 onto mesh x1
1727 // Note: u1 is defined so that mass conserved locally in each x1-elem
1728 //
1729 // Cao: 07/10/2009
1730 //
1732 
1733 void US_LammAstfvm::fun_phi( double x, double* y )
1734 {
1735  double x2 = x * x;
1736  y[ 0 ] = 0.5 * ( x2 - x );
1737  y[ 1 ] = 1.0 - x2;
1738  y[ 2 ] = 0.5 * ( x2 + x );
1739 }
1740 
1741 void US_LammAstfvm::fun_dphi( double x, double* y )
1742 {
1743  // quadratic basis
1744  double x2 = x * x;
1745  y[ 0 ] = 0.5 * ( x2 - x );
1746  y[ 1 ] = 1.0 - x2;
1747  y[ 2 ] = 0.5 * ( x2 + x );
1748 
1749  // derivatives
1750  y[ 3 ] = x - 0.5;
1751  y[ 4 ] = -2. * x;
1752  y[ 5 ] = x + 0.5;
1753 }
1754 
1755 void US_LammAstfvm::fun_Iphi( double x, double* y )
1756 {
1757  double x2 = x * x;
1758  double x3 = x * x2 / 6.0;
1759  x2 *= 0.25;
1760  y[ 0 ] = x3 - x2;
1761  y[ 1 ] = x - x3 * 2.0;
1762  y[ 2 ] = x3 + x2;
1763 }
1764 
1766 //
1767 // integrate a quadratic function defined on (x[0], x[1])
1768 // by nodal values u[0], u[1], u[2] from xi=xia to xib
1769 // here x=x[0]+(x[1]-x[0])*(xi+1)/2
1770 //
1772 double US_LammAstfvm::IntQ( double *x, double *u, double xia, double xib )
1773 {
1774  double intgrl;
1775  double phia[ 3 ];
1776  double phib[ 3 ];
1777 
1778  fun_Iphi( xia, phia );
1779  fun_Iphi( xib, phib );
1780 
1781  intgrl = ( x[ 1 ] - x[ 0 ] ) / 2. * (
1782  u[ 0 ] * ( phib[ 0 ] - phia[ 0 ] ) +
1783  u[ 1 ] * ( phib[ 1 ] - phia[ 1 ] ) +
1784  u[ 2 ] * ( phib[ 2 ] - phia[ 2 ] ) );
1785 
1786  return( intgrl );
1787 }
1788 
1790 //
1791 // integrate a piecewise quadratic function defined on mesh *x
1792 // by nodal values *u from xia in elem ka to xib in elem kb
1793 //
1795 double US_LammAstfvm::IntQs( double *x, double *u,
1796  int ka, double xia, int kb, double xib )
1797 {
1798  double intgrl;
1799 
1800  if ( ka == kb )
1801  {
1802  intgrl = IntQ( x + ka, u + 2 * ka, xia, xib );
1803  }
1804 
1805  else // integral across different elems of mesh x
1806  {
1807  intgrl = IntQ( x + ka, u + 2 * ka, xia, 1. );
1808  intgrl += IntQ( x + kb, u + 2 * kb, -1., xib );
1809 
1810  for ( int k = ka + 1; k <= kb - 1; k++ )
1811  {
1812  int k2 = k + k;
1813  intgrl += ( x[ k + 1 ] - x[ k ] )
1814  * ( u[ k2 ] + u[ k2 + 2 ] + 4. * u[ k2 + 1 ] ) / 6.;
1815  }
1816  }
1817 
1818  return ( intgrl );
1819 }
1820 
1822 //
1823 // Interpolation-projection of a piecewise quadratic u0 on x0
1824 // onto piecewise quadratic on mesh x1
1825 //
1826 // M0, M1 = number of elems in x0 and x1
1827 //
1829 void US_LammAstfvm::ProjectQ( int M0, double *x0, double *u0,
1830  int M1, double *x1, double *u1 )
1831 {
1832  double intgrl;
1833  double phi[ 3 ];
1834 
1835  int* ke = new int [ M1 + 1 ];
1836  double* xi = new double [ M1 + 1 ];
1837 
1838  LocateStar( M0 + 1, x0, M1 + 1, x1, ke, xi );
1839 
1840  // u1 = u0 at all nodes
1841  for ( int j = 0; j <= M1; j++ )
1842  {
1843  fun_phi( xi[ j ], phi );
1844 
1845  int idx = 2 * ke[ j ];
1846  u1[ 2 * j ] = phi[ 0 ] * u0[ idx ]
1847  + phi[ 1 ] * u0[ idx + 1 ]
1848  + phi[ 2 ] * u0[ idx + 2 ];
1849  }
1850 
1851  for ( int j = 0; j < M1; j++ )
1852  {
1853  int j2 = 2 * j;
1854 
1855  intgrl = IntQs( x0, u0, ke[ j ], xi[ j ], ke[ j + 1 ], xi[ j + 1 ] );
1856 
1857  u1[ j2 + 1 ] = 1.5 * intgrl / ( x1[ j + 1 ] - x1[ j ] )
1858  - 0.25 * ( u1[ j2 ] + u1[ j2 + 2 ] );
1859  }
1860 
1861  delete [] ke;
1862  delete [] xi;
1863 }
1864 
1865 
1867 //
1868 // LsSolver53
1869 // Matrix A[m+1][5]: with A[i][2]=a_ii diagonal
1870 // A[i][0]=A[i][4]=0, for i = odd
1871 // m must be even
1872 //
1874 void US_LammAstfvm::LsSolver53( int m, double **A, double *b, double *x )
1875 {
1876 
1877  for ( int j = 0; j < m - 1; j += 2 )
1878  {
1879  int j1 = j + 1;
1880  int j2 = j + 2;
1881 
1882  double multi = -A[ j1 ][ 1 ] / A[ j ][ 2 ];
1883  A[ j1 ][ 2 ] += multi * A[ j ][ 3 ];
1884  A[ j1 ][ 3 ] += multi * A[ j ][ 4 ];
1885  b[ j1 ] += multi * b[ j ];
1886 
1887  multi = -A[ j2 ][ 0 ] / A[ j ][ 2 ];
1888  A[ j2 ][ 1 ] += multi * A[ j ][ 3 ];
1889  A[ j2 ][ 2 ] += multi * A[ j ][ 4 ];
1890  b[ j2 ] += multi * b[ j ];
1891 
1892  multi = -A[ j2 ][ 1 ] / A[ j1 ][ 2 ];
1893  A[ j2 ][ 2 ] += multi * A[ j1 ][ 3 ];
1894  b[ j2 ] += multi * b[ j1 ];
1895  }
1896 
1897  // Back-substitution
1898  x[ m ] = b[ m ] / A[ m ][ 2 ];
1899 
1900  for ( int j = m - 1; j > 0; j -= 2 )
1901  {
1902  int jm = j - 1;
1903  int jp = j + 1;
1904 
1905  x[ j ] = ( b[ j ] - A[ j ][ 3 ] * x[ jp ] )
1906  / A[ j ][ 2 ];
1907 
1908  x[ jm ] = ( b[ jm ] - A[ jm ][ 3 ] * x[ j ]
1909  - A[ jm ][ 4 ] * x[ jp ] )
1910  / A[ jm ][ 2 ];
1911  }
1912 
1913 }
1914 
1915 // determine the non-ideal case number: 0/1/2/3
1916 int US_LammAstfvm::nonIdealCaseNo()
1917 {
1918  int rc = 0;
1919 
1920  if ( comp_x > 0 )
1921  return rc;
1922 
1923  NonIdealCaseNo = 0; // ideal
1924 
1925  US_Model::SimulationComponent* sc = &model.components[ comp_x ];
1926 
1927  if ( sc->sigma != 0.0 || sc->delta != 0.0 )
1928  { // non-zero sigma or delta given: concentration-dependent
1929  NonIdealCaseNo = 1;
1930  }
1931 
1932  if ( model.coSedSolute >= 0 )
1933  { // co-sedimentation solute index not -1: co-sedimenting
1934  if ( NonIdealCaseNo != 0 )
1935  rc = 2;
1936 
1937  NonIdealCaseNo = 2;
1938  }
1939 
1940  if ( compressib > 0.0 )
1941  { // compressibility factor positive: compressibility
1942  if ( NonIdealCaseNo != 0 )
1943  rc = 3;
1944 
1945  NonIdealCaseNo = 3;
1946  }
1947  return rc;
1948 }
1949 
1950 // perform quadratic interpolation to fill out full concentration vector
1951 void US_LammAstfvm::quadInterpolate( double* x0, double* u0, int N0,
1952  QVector< double >& xout, QVector< double >& cout )
1953 {
1954  int nout = xout.size(); // output concentrations count
1955  int kk = 0; // initial output index
1956  double xv = xout[ 0 ]; // first output X
1957  double yv; // output Y
1958 
1959  int ii = 2; // next x0 index to x3
1960  int jj = 3; // next u0 index to y2
1961  double x1 = x0[ 0 ]; // initial start X
1962  double x3 = x0[ 1 ]; // initial end X
1963  double x2 = ( x1 + x3 ) * 0.5; // initial mid-point X
1964  double y1 = u0[ 0 ]; // initial start Y
1965  double y2 = u0[ 1 ]; // initial mid-point Y
1966  double y3 = u0[ 2 ]; // initial end Y
1967 
1968  cout.resize( nout );
1969 
1970  while ( kk < nout )
1971  { // loop to output interpolated concentrations
1972  xv = xout[ kk ]; // X for which we need a Y
1973 
1974  while ( xv > x3 && ii < N0 )
1975  { // if need be, walk up input until between x values
1976  x1 = x3; // start x (previous end)
1977  x3 = x0[ ii++ ]; // end (next) x
1978  y1 = y3; // y at start x
1979  y2 = u0[ jj++ ]; // y at mid-point
1980  y3 = u0[ jj++ ]; // y at end (next) x
1981 
1982  }
1983  x2 = ( x1 + x3 ) * 0.5; // mid-point x
1984 //y1 /= x1;
1985 //y2 /= x2;
1986 //y3 /= x3;
1987 
1988  // do the quadratic interpolation of this Y (C*r)
1989  yv =
1990  (( ( xv - x2 ) * ( xv - x3 ) ) / ( ( x1 - x2 ) * ( x1 - x3 ) )) * y1 +
1991  (( ( xv - x1 ) * ( xv - x3 ) ) / ( ( x2 - x1 ) * ( x2 - x3 ) )) * y2 +
1992  (( ( xv - x1 ) * ( xv - x2 ) ) / ( ( x3 - x1 ) * ( x3 - x2 ) )) * y3;
1993 //y1 *= x1;
1994 //y2 *= x2;
1995 //y3 *= x3;
1996 
1997  // output interpolated concentration with r factor removed (C = (C*r)/r)
1998  cout[ kk++ ] = yv / xv;
1999 //cout[ kk++ ] = yv;
2000  }
2001 
2002 }
2003 
2004 // load MfemData object used internally from caller's RawData object
2005 void US_LammAstfvm::load_mfem_data( US_DataIO::RawData& edata,
2006  US_AstfemMath::MfemData& fdata,
2007  bool zeroout )
2008 {
2009  int nscan = edata.scanData.size(); // scan count
2010 // int nconc = edata.x.size(); // concentrations count
2011  int nconc = edata.xvalues.size(); // concentrations count
2012  fdata.id = edata.description;
2013  fdata.cell = edata.cell;
2014  fdata.scan .resize( nscan ); // mirror number of scans
2015  fdata.radius.resize( nconc ); // mirror number of radius values
2016 
2017  for ( int ii = 0; ii < nscan; ii++ )
2018  { // copy over all scans
2019  US_AstfemMath::MfemScan* fscan = &fdata.scan[ ii ];
2020 
2021  fscan->temperature = edata.scanData[ ii ].temperature;
2022  fscan->rpm = edata.scanData[ ii ].rpm;
2023  fscan->time = edata.scanData[ ii ].seconds;
2024  fscan->omega_s_t = edata.scanData[ ii ].omega2t;
2025 #if 0
2026  fscan->conc.resize( nconc ); // mirror number of concentrations
2027 
2028  if ( zeroout )
2029  { // if so specified, set concentrations to zero
2030  fscan->conc.fill( 0.0 );
2031  }
2032 
2033  else
2034  { // otherwise, copy input
2035  for ( int jj = 0; jj < nconc; jj++ )
2036  { // copy all concentrations for a scan
2037  fscan->conc[ ii ] = edata.value( ii, jj );
2038  }
2039  }
2040  }
2041 
2042  for ( int jj = 0; jj < nconc; jj++ )
2043  { // copy all radius values
2044  fdata.radius[ jj ] = edata.radius( jj );
2045  }
2046 #endif
2047 #if 1
2048  if ( zeroout )
2049  fscan->conc.fill( 0.0, nconc ); // if so specified, set concentrations to zero
2050 
2051  else // Otherwise, copy concentrations
2052  fscan->conc = edata.scanData[ ii ].rvalues;
2053  }
2054  fdata.radius = edata.xvalues;
2055 int nn=fdata.radius.size() - 1;
2056 int mm=nn/2;
2057 DbgLv(2) << "LdDa: n r0 rm rn" << nn << fdata.radius[0] << fdata.radius[mm] << fdata.radius[nn];
2058 #endif
2059 }
2060 
2061 // store MfemData object used internally into caller's RawData object
2062 void US_LammAstfvm::store_mfem_data( US_DataIO::RawData& edata,
2063  US_AstfemMath::MfemData& fdata )
2064 {
2065  int nscan = fdata.scan.size(); // scan count
2066  int nconc = fdata.radius.size(); // concentrations count
2067 
2068  edata.description = fdata.id;
2069  edata.cell = fdata.cell;
2070  edata.xvalues = fdata.radius;
2071  edata.scanData.resize( nscan ); // mirror number of scans
2072 
2073  for ( int ii = 0; ii < nscan; ii++ )
2074  { // copy over each scan
2075  US_AstfemMath::MfemScan* fscan = &fdata.scan [ ii ];
2076  US_DataIO::Scan* escan = &edata.scanData[ ii ];
2077 
2078  escan->temperature = fscan->temperature;
2079  escan->rpm = fscan->rpm;
2080  escan->seconds = fscan->time;
2081  escan->omega2t = fscan->omega_s_t;
2082  escan->plateau = fdata.radius[ nconc - 1 ];
2083  escan->rvalues = fscan->conc;
2084  }
2085 }
2086 
2087 void US_LammAstfvm::setStopFlag( bool flag )
2088 {
2089  stopFlag = flag;
2090  qApp->processEvents();
2091 DbgLv(2) << "setStopFlag" << stopFlag;
2092 }
2093 
2094 void US_LammAstfvm::setMovieFlag( bool flag )
2095 {
2096  movieFlag = flag;
2097  qApp->processEvents();
2098 DbgLv(2) << "setMovieFlag" << movieFlag;
2099 }
2100 
Generated on Sun Jun 21 2015 11:40:42 for UltraScan III by   doxygen 1.8.3.1-20130324