UltraScan III
us_astfem_rsa.cpp
Go to the documentation of this file.
1 
3 #include "us_astfem_rsa.h"
4 #include "us_astfem_math.h"
5 #include "us_hardware.h"
6 #include "us_math2.h"
7 #include "us_memory.h"
8 #include "us_stiffbase.h"
9 #include "us_settings.h"
10 #include "us_sleep.h"
11 #include "us_util.h"
12 
13 US_Astfem_RSA::US_Astfem_RSA( US_Model& model,
14  US_SimulationParameters& params,
15  QObject* parent )
16  : QObject( parent ), system( model ), simparams( params )
17 {
18  stopFlag = false;
19  use_time = false;
20  time_correction = true;
21  simout_flag = false;
22  show_movie = false;
23  dbg_level = 0;
24 }
25 
26 int US_Astfem_RSA::calculate( US_DataIO::RawData& exp_data )
27 {
28  US_AstfemMath::MfemInitial* vC0 = NULL; // Initial concentration for
29  // multiple components
30  US_AstfemMath::MfemInitial CT0; // Initial total concentration
31  US_AstfemMath::MfemData simdata;
32  double current_time = 0.0;
33  double current_om2t = 0.0;
34  double current_speed = 0.0;
35  double duration = 0.0;
36  double delay = 0.0;
37 
38  int initial_npts = 1000;
39  int current_assoc = 0;
40  int size_cv = system.components.size();
41  QVector< bool > reactVec( size_cv );
42  bool* reacting = reactVec.data();
43 
44  double accel_time;
45  double dr;
46 #if 0
47 #define TIMING_RA_INC 500
48 #define TIMING_RA
49 #endif
50 #ifdef TIMING_RA
51 QDateTime calcStart = QDateTime::currentDateTime();
52 static int ncalls=0;
53 static int totTC=0;
54 static int totT1=0;
55 static int totT2=0;
56 static int totT3=0;
57 static int totT4=0;
58 static int totT5=0;
59 static int totT6=0;
60 static int totT7=0;
61 static int totT8=0;
62 #endif
63 
64  af_params.simpoints = simparams.simpoints;
65  af_params.pathlength = simparams.cp_pathlen;
66  af_params.dt = 1.0;
67  af_params.time_steps = simparams.simpoints;
68  af_params.first_speed = simparams.speed_step[ 0 ].rotorspeed;
69  af_params.omega_s = sq( (double)af_params.first_speed * M_PI / 30.0 );
70  af_params.start_om2t = af_params.omega_s;
71  //af_params.start_time = 0.0;
72  af_params.current_meniscus = simparams.meniscus;
73  af_params.current_bottom = simparams.bottom;
74 
75  load_mfem_data( exp_data, af_data );
76 
77  int npts = af_data.scan[ 0 ].conc.size();
78  initial_npts = ( initial_npts < npts ) ? initial_npts : npts;
79 DbgLv(2) << "RSA: AP.rotorspeed" << af_params.first_speed;
80 DbgLv(2) << "RSA: AP.simpoints" << af_params.simpoints;
81 DbgLv(2) << "RSA: AP.pathlength" << af_params.pathlength;
82 DbgLv(2) << "RSA: AP.time_steps" << af_params.time_steps;
83 DbgLv(2) << "RSA: AP.first_speed" << af_params.first_speed;
84 DbgLv(2) << "RSA: AP.omega_s" << af_params.omega_s;
85 DbgLv(2) << "RSA: AP.current_meniscus" << af_params.current_meniscus;
86 DbgLv(2) << "RSA: AP.current_bottom" << af_params.current_bottom;
87 DbgLv(2) << "RSA: AP.start_om2t" << af_params.start_om2t;
88 DbgLv(2) << "RSA: scan0size" << npts;
89 DbgLv(2) << "RSA: af_c0size" << initial_npts;
90 // update_assocv();
91  initialize_rg(); // Reaction group
92 //*DEBUG*
93 if(dbg_level>1)
94 {
95  int nrg=rg.size();
96  DbgLv(2) << "RSA: RG size" << nrg;
97  for(int jr=0;jr<nrg;jr++)
98  {
99  int na=rg[jr].association.size();
100  int ng=rg[jr].GroupComponent.size();
101  DbgLv(2) << "RSA: RG" << jr << ": na" << na << "ng" << ng;
102  for(int ja=0;ja<na;ja++)
103  DbgLv(2) << "RSA: Assoc" << ja << ":" << rg[jr].association[ja];
104  for(int jg=0;jg<ng;jg++)
105  DbgLv(2) << "RSA: GrpComp" << jg << ":" << rg[jr].GroupComponent[jg];
106  }
107 }
108 //*DEBUG*
109  adjust_limits( af_params.first_speed );
110 DbgLv(2) << "RSA: sbottom acbottom" << simparams.bottom
111  << af_params.current_bottom;
112  af_data.meniscus = af_params.current_meniscus;
113  af_data.bottom = af_params.current_bottom;
114 // cv = component vector
115  for ( int cc = 0; cc < size_cv; cc++ )
116  {
117 #ifdef TIMING_RA
118 QDateTime clcSt1 = QDateTime::currentDateTime();
119 #endif
120  US_Model::SimulationComponent* sc = &system.components[ cc ];
121  US_Model::Association* as;
122  reacting[ cc ] = false;
123 DbgLv(2) << "RSA: cc assoc.size" << cc << system.associations.size();
124 
125  for ( int aa = 0; aa < system.associations.size(); aa++ )
126  {
127  as = &system.associations[ aa ];
128 DbgLv(2) << "RSA: aa react.size" << aa << as->rcomps.size();
129 
130  for ( int jj = 0; jj < as->rcomps.size(); jj++ )
131  {
132  if ( cc == (int)as->rcomps[ jj ] )
133  {
134  reacting[ cc ] = true;
135  current_assoc = aa;
136  break; // Since a comp appears at most once in an assoc rule
137  }
138  }
139 DbgLv(2) << "RSA: cc aa current_assoc" << cc << aa << current_assoc;
140  }
141 
142  current_time = 0.0;
143  current_om2t = 0.0;
144  last_time = 0.0;
145  current_speed = 0.0;
146  w2t_integral = 0.0;
147 
148  CT0.radius .clear();
149  CT0.concentration.clear();
150  CT0.radius .reserve( initial_npts );
151  CT0.concentration.reserve( initial_npts );
152 
153  dr = ( af_params.current_bottom - af_params.current_meniscus ) /
154  ( initial_npts - 1 );
155 
156  for ( int jj = 0; jj < initial_npts; jj++ )
157  {
158  CT0.radius.append( af_params.current_meniscus + jj * dr );
159  CT0.concentration.append( 0.0 );
160  //CT0.concentration.append( sc->signal_concentration );
161  }
162 
163  af_c0.radius .clear();
164  af_c0.concentration.clear();
165 
166  // Once time invariant noise has been removed in a band experiment, we
167  // can use the first scan of the first speed step of the experiment as
168  // the initial concentration of the simulation. The approach will copy
169  // the 1st scan concentration vector into each component's c0 vector.
170  // NNLS will scale the appropriate concentration for each component. The
171  // assumption is made that any potential differentiation of components in
172  // the initial scan is minimal compared to any solute flow disturbances
173  // at the meniscus. For this approach to work well it is necessary to
174  // pick the first data point close to the meniscus and to include the
175  // earliest possible scan in the experiment. Also, time invariant noise
176  // should be subtracted first.
177 
178  if ( simparams.firstScanIsConcentration )
179  {
180  US_AstfemMath::MfemScan* scan0 = &af_data.scan[ 0 ];
181  US_AstfemMath::MfemInitial scan1;
182  scan1.radius .clear();
183  scan1.concentration.clear();
184  scan1.radius .reserve( af_data.radius.size() );
185  scan1.concentration.reserve( af_data.radius.size() );
186 
187  for ( int jj = 0; jj < af_data.radius.size(); jj++ )
188  {
189  scan1.radius .append( af_data.radius[ jj ] );
190  scan1.concentration.append( scan0->conc [ jj ] );
191  }
192 
193  US_AstfemMath::interpolate_C0( scan1, af_c0 );
194  }
195 
196  if ( ! reacting[ cc ] ) // noninteracting
197  {
198  initialize_conc( cc, CT0, true );
199 int nc0=CT0.concentration.size()-1;
200 int mc0=nc0/2;
201 DbgLv(2) << "RSA: s0 i_conc cc step lsc" << cc << 0 << 0 << ": c0 cm cn"
202  << CT0.concentration[0] << CT0.concentration[mc0] << CT0.concentration[nc0];
203 
204  af_params.s .resize( 1 );
205  af_params.D .resize( 1 );
206  af_params.kext.resize( 1 );
207 
208  af_params.s [ 0 ] = sc->s;
209  af_params.D [ 0 ] = sc->D;
210  //af_params.kext[ 0 ] = sc->extinction;
211  af_params.kext[ 0 ] = sc->extinction * af_params.pathlength;
212 #ifdef TIMING_RA
213 totT1+=(clcSt1.msecsTo(QDateTime::currentDateTime()));
214 #endif
215  int lscan = 0;
216  int fscan = 0;
217  int nstep = simparams.speed_step.size();
218  double time0 = 0.0;
219  double time1 = 0.0;
220  double time2 = 0.0;
221  double omeg0 = 0.0;
222  double omeg1 = 0.0;
223  double omeg2 = 0.0;
224  double step_speed;
225  US_SimulationParameters::SpeedProfile* sp;
226  US_AstfemMath::MfemData* ed;
227 
228  for ( int step = 0; step < nstep; step++ )
229  {
230 #ifdef TIMING_RA
231 QDateTime clcSt2 = QDateTime::currentDateTime();
232 #endif
233  sp = &simparams.speed_step[ step ];
234  ed = &af_data;
235  step_speed = (double)sp->rotorspeed;
236  fscan = 0;
237  lscan = fscan + af_data.scan.count() - 1;
238  bool in_step = ( ed->scan[ fscan ].time <= sp->time_last &&
239  ed->scan[ lscan ].time >= sp->time_first );
240 
241 DbgLv(2) << "RSA:PO: STEP NSTEP" << step << nstep
242  << "ETIMEF ETIMEL" << ed->scan[fscan].time << ed->scan[lscan].time
243  << "STIMEF STIMEL" << sp->time_first << sp->time_last;
244  // For multi-speed, break out once speed step is beyond data
245  if ( step > 0 && ed->scan[ lscan ].time < sp->time_first )
246  break;
247 
248  adjust_limits( sp->rotorspeed );
249 
250  time0 = time2;
251  omeg0 = omeg2;
252  time1 = sp->time_first;
253  time2 = sp->time_last;
254  omeg1 = sp->w2t_first;
255  omeg2 = sp->w2t_last;
256  ed->meniscus = af_params.current_meniscus;
257  ed->bottom = af_params.current_bottom;
258  accel_time = 0.0;
259 DbgLv(2) << "RSA:PO: FSCAN TIME1" << fscan << time1 << "cc step" << cc << step
260  << "menisc bott" << ed->meniscus << ed->bottom;
261 
262  // We need to simulate acceleration
263  if ( sp->acceleration_flag )
264  {
265  // If the speed difference is larger than acceleration rate then
266  // we have at least 1 acceleration step
267 
268  af_params.time_steps = qRound(
269  qAbs( step_speed - current_speed ) / sp->acceleration );
270 
271  // Each simulation step is 1 second in the acceleration phase
272  // Use a fixed grid with refinement at both ends and with
273  // twice the number of points
274  af_params.dt = 1.0;
275  af_params.simpoints = 2 * simparams.simpoints;
276 DbgLv(2) << "RSA: AP during ACCELERATION";
277 DbgLv(2) << "RSA: AP.rotorspeed" << af_params.first_speed;
278 DbgLv(2) << "RSA: AP.simpoints" << af_params.simpoints;
279 DbgLv(2) << "RSA: AP.pathlength" << af_params.pathlength;
280 DbgLv(2) << "RSA: AP.time_steps" << af_params.time_steps;
281 DbgLv(2) << "RSA: AP.first_speed" << af_params.first_speed;
282 DbgLv(2) << "RSA: AP.omega_s" << af_params.omega_s;
283 DbgLv(2) << "RSA: AP.current_meniscus" << af_params.current_meniscus;
284 DbgLv(2) << "RSA: AP.current_bottom" << af_params.current_bottom;
285 DbgLv(2) << "RSA: AP.start_om2t" << af_params.start_om2t;
286 DbgLv(2) << "RSA: scan0size" << npts;
287 DbgLv(2) << "RSA: af_c0size" << initial_npts;
288  accel_time = af_params.dt * af_params.time_steps;
289  delay = time1 - time0;
290 
291  // Do calculations to find the position of the acceleration
292  // zone within the gap between speed steps.
293  //
294  // time0 t2 t3 time1
295  // | sp1 | sp1->sp2 | sp2 |
296  // | dt1 | dt2 | dt3 |
297  // -------------------------------
298  // | dw1 | dw2 | dw3 |
299  // omeg0 w2 w3 omeg1
300  //
301  // The code below is meant to determine "dt1", the time
302  // difference between the end of the previous step and the
303  // beginning of the acceleration zone. This is done by
304  // solving the simultaneous equations:
305  // dt1 * ddw1 + dt3 * ddw3 = omeg1 - omeg0 - dw2
306  // dt1 + dt3 = time1 - time0 - dt2
307  // which leads to:
308  //
309  // dt1 = ( omeg1 - omeg0 - dw2 - drtm * ddw3 )
310  // / ( ddw1 - ddw3 )
311  //
312  // for ddw1, ddw3 being changes in omega per second in 2 zones
313  // and drtm = time1 - time0 - dt2
314 
315  double drtm = delay - accel_time;
316  double wfac = af_params.dt * sq( M_PI / 30.0 );
317 DbgLv(2) << "RSA:PO: time0 time1 omeg0 omeg1"
318  << time0 << time1 << omeg0 << omeg1 << "delay atime" << delay << accel_time;
319  double dw2 = 0.0;
320  double dt2 = accel_time / af_params.dt;
321  double ddw1 = wfac * sq( current_speed );
322  double ddw3 = wfac * sq( step_speed );
323  int ndt2 = qRound( dt2 );
324  double crpm = current_speed;
325  double rpmi = ( step_speed - current_speed ) / dt2;
326 
327  for ( int kk = 0; kk < ndt2; kk++ )
328  { // Accumulate omega^2t over acceleration zone
329  crpm += rpmi;
330  dw2 += ( wfac * sq( crpm ) );
331  }
332 
333  double dw1 = 0.0;
334  double dt1 = 0.0;
335 
336  if ( current_speed == 0.0 )
337  { // For the first (only?) speed step
338  double dt3 = ( omeg1 - dw2 ) / ddw3;
339  dt1 = (double)qFloor( time1 - dt3 - accel_time );
340  dt1 = ( dt1 < 0.0 ) ? 1.0 : dt1;
341  dw1 = wfac * sq( rpmi );
342  }
343 
344  else
345  { // For speed steps beyond the first
346  dt1 = (double)qFloor(
347  ( omeg1 - omeg0 - dw2 - drtm * ddw3 )
348  / ( ddw1 - ddw3 ) );
349  dw1 = ddw1 * dt1;
350  }
351 DbgLv(2) << "RSA:PO: ddw1 ddw3" << ddw1 << ddw3 << "dt2 dw2" << dt2 << dw2
352  << "dt1 dw1" << dt1 << dw1;
353 
354  af_params.start_time = time0 + dt1;
355  af_params.start_om2t = omeg0 + dw1;
356  int nradi = simdata.radius.size();
357 DbgLv(2) << "RSA:PO: st_ time om2t"
358  << af_params.start_time << af_params.start_om2t;
359 nc0=CT0.concentration.size()-1;
360 mc0=nc0/2;
361 DbgLv(2) << "RSA: S:cc step" << cc << step << ": c0 cm cn"
362  << CT0.concentration[0] << CT0.concentration[mc0] << CT0.concentration[nc0];
363 
364  // If beyond the 1st speed step and acceleration begins
365  // a little after the end of the previous step, calculate
366  // the simulation for constant speed leading up to the
367  // acceleration zone
368  if ( step > 0 && dt1 > 2.0 )
369 // if ( step > 0 && dt1 > 1.0 )
370  {
371  int svnpts = af_params.time_steps;
372  af_params.start_time = time0;
373  af_params.start_om2t = omeg0;
374 // af_params.time_steps = qRound( dt1 ) + 1;
375  af_params.time_steps = qRound( dt1 );
376  dt1 = (double)af_params.time_steps;
377  dw1 = ddw1 * dt1;
378 
379  adjust_limits( qRound( current_speed ) );
380 
381  calculate_ni( current_speed, current_speed,
382  CT0, simdata, false );
383 
384  af_params.start_time += dt1;
385  af_params.start_om2t += dw1;
386  af_params.time_steps = svnpts;
387 DbgLv(2) << "RSA:PO: Accel nradi" << nradi << "CT0size" << CT0.radius.size();
388  }
389 
390  // Calculate the simulation for the acceleration zone
391 
392 // adjust_limits( qRound( ( current_speed + step_speed ) * 0.5 ) );
393  adjust_limits( sp->rotorspeed );
394 
395  calculate_ni( current_speed, step_speed,
396  CT0, simdata, true );
397 
398  qApp->processEvents();
399  if ( stopFlag ) return 1;
400 
401  // Add the acceleration time:
402  current_time = af_params.start_time + accel_time;
403  current_om2t = af_params.start_om2t + dw2;
404 
405 #ifndef NO_DB
406  emit new_time( current_time );
407  qApp->processEvents();
408 #endif
409 // adjust_limits( sp->rotorspeed );
410  } // End of acceleration
411 
412  else
413  { // No acceleration
414  current_time = time1;
415  current_om2t = omeg1;
416  }
417 #ifdef TIMING_RA
418 QDateTime clcSt3 = QDateTime::currentDateTime();
419 totT2+=(clcSt2.msecsTo(clcSt3));
420 #endif
421 
422  duration = time2 - current_time;
423 
424  if ( step == simparams.speed_step.size() - 1 )
425  duration += (double)( (int)( duration * 0.05 ) ); // +5%
426 
427  if ( accel_time > duration )
428  {
429  DbgErr() << "Attention: acceleration time exceeds duration - "
430  "please check initialization\n";
431  return -1;
432  }
433 
434  double omega = step_speed * M_PI / 30.0;
435  af_params.omega_s = sq( omega );
436 DbgLv(2) << "RSA: ***COMPUTED omega_s:" << af_params.omega_s << "omega step_speed"
437  << omega << step_speed;
438 
439  double lg_bm_rat = log(
440  af_params.current_bottom / af_params.current_meniscus );
441  double s_omg_fac = qAbs( sc->s ) * af_params.omega_s;
442  af_params.dt = lg_bm_rat
443  / ( s_omg_fac * ( simparams.simpoints - 1 ) );
444 
445  if ( af_params.dt > duration)
446  {
447 double dtsv=af_params.dt;
448  af_params.dt = duration;
449  af_params.simpoints = 1 +
450  (int)( lg_bm_rat / ( s_omg_fac * af_params.dt ) );
451 DbgLv(1) << "RSA: ***CORRECTED dt:" << duration << dtsv << af_params.simpoints;
452 DbgLv(1) << "RSA: ***CORRECTED dt: lg_bm_rat" << lg_bm_rat << "sc->s" << sc->s
453  << "omega_s" << af_params.omega_s << "omg_fac" << s_omg_fac;
454  }
455 
456  if ( af_params.simpoints > 10000 )
457  {
458 DbgLv(1) << "RSA: ***COMPUTED simpoints:" << af_params.simpoints
459  << "omg2t s_val dt dura" << af_params.omega_s << sc->s << af_params.dt
460  << duration << "** simpts set to 10000 **";
461  af_params.simpoints = 10000;
462  }
463 
464  // Find out the minimum number of simpoints needed to provide
465  // the necessary dt:
466 // af_params.time_steps = (int)( duration / af_params.dt ) + 1;
467 // af_params.time_steps = qRound( duration / af_params.dt ) + 1;
468  af_params.time_steps = qRound( duration / af_params.dt );
469  af_params.start_time = current_time;
470  af_params.start_om2t = current_om2t;
471 DbgLv(2) << "RSA: (2)AP.start_om2t" << af_params.start_om2t;
472 // af_params.dt = duration / ( af_params.time_steps - 1 );
473  af_params.dt = duration / (double)af_params.time_steps;
474 if(af_params.time_steps>10000)
475 DbgLv(1) << "RSA: ***COMPUTED time_steps:" << af_params.time_steps
476  << "dur" << duration << "o2t s" << af_params.omega_s << sc->s
477  << "dt" << af_params.dt << "rat fac" << lg_bm_rat << s_omg_fac;
478 #ifdef TIMING_RA
479 QDateTime clcSt4 = QDateTime::currentDateTime();
480 totT3+=(clcSt3.msecsTo(clcSt4));
481 #endif
482 // current_time = time1;
483 nc0=CT0.concentration.size()-1;
484 mc0=nc0/2;
485 DbgLv(2) << "RSA: B:cc step" << cc << step << ": c0 cm cn"
486  << CT0.concentration[0] << CT0.concentration[mc0] << CT0.concentration[nc0];
487 US_AstfemMath::MfemData* ed = &af_data;
488 int nr1=ed->radius.size();
489 DbgLv(2) << "RSA:(1) nr1" << nr1 << "r sme"
490  << ed->radius[0] << ed->radius[1] << ed->radius[2]
491  << ed->radius[nr1/2-1] << ed->radius[nr1/2] << ed->radius[nr1/2+1]
492  << ed->radius[nr1-3] << ed->radius[nr1-2] << ed->radius[nr1-1];
493 
494  // Calculate the simulation for the bulk of the speed step
495 
496  calculate_ni( step_speed, step_speed,
497  CT0, simdata, false );
498 
499  qApp->processEvents();
500  if ( stopFlag ) return 1;
501 
502  // Set the current time to the last scan of this speed step
503 // duration = time2 - time0;
504 // duration = time2 - current_time;
505  current_time = time2;
506  current_om2t = omeg2;
507 DbgLv(2) << "RSA:T step rpm" << step << sp->rotorspeed
508  << "st_ curr_time" << af_params.start_time << current_time
509  << "fscan lscan" << fscan << lscan << "dt" << af_params.dt;
510 int mmm=simdata.scan.size()-1;
511 int kkk=qMin(af_params.time_steps-1,mmm);
512 DbgLv(2) << "RSA:T eomg1 eomg2" << ed->scan[fscan].omega_s_t
513  << ed->scan[lscan].omega_s_t << " somg1 somg2"
514  << simdata.scan[0].omega_s_t << simdata.scan[mmm].omega_s_t;
515 DbgLv(2) << "RSA:T etim1 etim2" << ed->scan[fscan].time
516  << ed->scan[lscan].time << " stim1 stim2"
517  << simdata.scan[0].time << simdata.scan[mmm].time << simdata.scan[kkk].time
518  << "sssz tsteps" << simdata.scan.size() << af_params.time_steps;
519 #ifdef TIMING_RA
520 QDateTime clcSt5 = QDateTime::currentDateTime();
521 totT4+=(clcSt4.msecsTo(clcSt5));
522 #endif
523 
524  // Interpolate the simulated data onto
525  // the experimental time and radius grid,
526  // if the experimental time range fits in this speed step.
527 int nr2=ed->radius.size();
528 DbgLv(2) << "RSA:(2) nr2" << nr2 << "r sme"
529  << ed->radius[0] << ed->radius[1] << ed->radius[2]
530  << ed->radius[nr2/2-1] << ed->radius[nr2/2] << ed->radius[nr2/2+1]
531  << ed->radius[nr2-3] << ed->radius[nr2-2] << ed->radius[nr2-1];
532  if ( in_step )
533  {
534 // US_AstfemMath::interpolate( *ed, simdata, use_time,
535 // fscan, ++lscan );
536  US_AstfemMath::interpolate( *ed, simdata, use_time );
537 DbgLv(2) << "RSA:T INTERP in step" << step;
538  }
539 int nr3=ed->radius.size();
540 DbgLv(2) << "RSA:(3) nr3" << nr3 << "r sme"
541  << ed->radius[0] << ed->radius[1] << ed->radius[2]
542  << ed->radius[nr3/2-1] << ed->radius[nr3/2] << ed->radius[nr3/2+1]
543  << ed->radius[nr3-3] << ed->radius[nr3-2] << ed->radius[nr3-1];
544 
545  if ( !simout_flag )
546  {
547  simdata.radius.clear();
548  simdata.scan .clear();
549  }
550 
551  // Set the current speed to the constant rotor speed of the
552  // current speed step
553  current_speed = step_speed;
554 
555 #ifndef NO_DB
556  qApp->processEvents();
557 #endif
558 #ifdef TIMING_RA
559 totT5+=(clcSt5.msecsTo(QDateTime::currentDateTime()));
560 #endif
561 DbgLv(2) << "RSA: step=" << step << "tsteps sttime" << af_params.time_steps
562  << current_time << "bottoms" << simparams.bottom_position
563  << simparams.bottom << af_params.current_bottom << "sp" << current_speed;
564 nc0=CT0.concentration.size()-1;
565 mc0=nc0/2;
566 DbgLv(2) << "RSA: E:cc step" << cc << step << ": c0 cm cn"
567  << CT0.concentration[0] << CT0.concentration[mc0] << CT0.concentration[nc0];
568 // af_params.first_speed = current_speed;
569 
570  } // Speed step loop
571 
572 #ifndef NO_DB
573  emit current_component( cc + 1 );
574  qApp->processEvents();
575 #endif
576  } // Non-interacting case
577  } // Model Component loop
578 #ifdef TIMING_RA
579 QDateTime clcSt6 = QDateTime::currentDateTime();
580 #endif
581 
582  // Resize af_params.local_index
583  af_params.local_index.resize( size_cv );
584 
585  US_AstfemMath::ComponentRole cr;
586 
587  for ( int group = 0; group < rg.size(); group++ )
588  {
589  int num_comp = rg[ group ].GroupComponent.size();
590  int num_rule = rg[ group ].association.size();
591  af_params.rg_index = group;
592 DbgLv(2) << "RSA: group nrule ncomp" << group << num_rule << num_comp;
593  af_params.s .resize( num_comp );
594  af_params.D .resize( num_comp );
595  af_params.kext .resize( num_comp );
596  af_params.role .resize( num_comp );
597  af_params.association.resize( num_rule );
598 
599  for ( int m = 0; m < num_rule; m++ )
600  {
601  af_params.association[ m ] =
602  system.associations[ rg[ group ].association[ m ] ];
603  }
604 
605  for ( int jj = 0; jj < num_comp; jj++ )
606  {
607  int index = rg[ group ].GroupComponent[ jj ];
608 DbgLv(2) << "RSA: jj index" << jj << index;
609 
610  US_Model::SimulationComponent* sc = &system.components[ index ];
611  af_params.s [ jj ] = sc->s;
612  af_params.D [ jj ] = sc->D;
613  af_params.kext[ jj ] = sc->extinction * af_params.pathlength;
614 //DbgLv(0) << "RSA: group jj" << group << jj << "extinct kext pathlen"
615 // << sc->extinction << af_params.kext[jj] << af_params.pathlength;
616 
617  // Global to local index
618  af_params.local_index[ index ] = jj;
619 
620  cr.comp_index = index;
621  cr.assocs .clear();
622  cr.stoichs.clear();
623 
624  af_params.role[ jj ] = cr; // Add jj'th rule
625 
626  // Check all assoc rule in this rg
627  for ( int aa = 0; aa < rg[ group ].association.size(); aa++ )
628  {
629  // Check all comp in rule
630  int rule = rg[ group ].association[ aa ];
631 DbgLv(2) << "RSA: aa rule" << aa << rule;
632  US_Model::Association* as = &system.associations[ rule ];
633 
634  for ( int rr = 0; rr < as->rcomps.size(); rr++ )
635  {
636 DbgLv(2) << "RSA: rr af-index as-react" << rr
637  << af_params.role[jj].comp_index << as->rcomps[rr];
638  if ( af_params.role[ jj ].comp_index ==
639  as->rcomps[ rr ] )
640  {
641  // local index for the rule
642  af_params.role[ jj ].assocs .append( aa );
643  af_params.role[ jj ].stoichs.append( as->stoichs[ rr ] );
644  break;
645  }
646  }
647  }
648  }
649 
650  for ( int aa = 0; aa < num_rule; aa++ )
651  {
652  US_Model::Association* as = &af_params.association[ aa ];
653  for ( int rr = 0; rr < as->rcomps.size(); rr++ )
654  {
655  as->rcomps[ rr ] =
656  af_params.local_index[ as->rcomps[ rr ] ];
657 DbgLv(2) << "RSA: aa rr rcn" << aa << rr << as->rcomps[rr];
658  }
659  }
660 
661  current_time = 0.0;
662  current_speed = 0.0;
663 DbgLv(2) << "RSA: (3)AP.start_om2t" << af_params.start_om2t;
664  w2t_integral = af_params.start_om2t;
665  last_time = af_params.start_time;
666 // w2t_integral = 0.0;
667 // last_time = 0.0;
668 
669  dr =
670  ( af_params.current_bottom - af_params.current_meniscus ) /
671  ( initial_npts - 1 );
672 
673  QVector< US_AstfemMath::MfemInitial > vC0Vec( num_comp );
674  vC0 = vC0Vec.data();
675 
676  for ( int jj = 0; jj < num_comp; jj++ )
677  {
678  CT0.radius .clear();
679  CT0.concentration.clear();
680  CT0.radius .reserve( initial_npts );
681  CT0.concentration.reserve( initial_npts );
682 DbgLv(2) << "RSA: jj in_npts" << jj << initial_npts;
683 
684  for ( int ii = 0; ii < initial_npts; ii++ )
685  {
686  CT0.radius .append( af_params.current_meniscus + ii * dr );
687  CT0.concentration.append( 0.0 );
688  }
689 
690  initialize_conc( rg[ group ].GroupComponent[ jj ], CT0, false );
691  vC0[ jj ] = CT0;
692  }
693 
694  decompose( vC0 );
695 DbgLv(2) << "RSA: decompose OUT";
696 #ifdef TIMING_RA
697 QDateTime clcSt5 = QDateTime::currentDateTime();
698 #endif
699  int nstep = simparams.speed_step.size();
700  int lscan = 0;
701  int fscan = 0;
702  double time0 = 0.0;
703  double time1 = 0.0;
704  double time2 = 0.0;
705  double step_speed;
706  US_SimulationParameters::SpeedProfile* sp;
707  US_AstfemMath::MfemData* ed;
708 ed=&af_data;
709 int nr4=ed->radius.size();
710 DbgLv(2) << "RSA:(4) nr4" << nr4 << "r sme"
711  << ed->radius[0] << ed->radius[1] << ed->radius[2]
712  << ed->radius[nr4/2-1] << ed->radius[nr4/2] << ed->radius[nr4/2+1]
713  << ed->radius[nr4-3] << ed->radius[nr4-2] << ed->radius[nr4-1];
714 
715  for ( int step = 0; step < nstep; step++ )
716  {
717  sp = &simparams.speed_step[ step ];
718  ed = &af_data;
719  step_speed = (double)sp->rotorspeed;
720 
721  adjust_limits( sp->rotorspeed );
722  ed->meniscus = af_params.current_meniscus;
723  ed->bottom = af_params.current_bottom;
724  accel_time = 0.0;
725  fscan = 0;
726  lscan = fscan + af_data.scan.count() - 1;
727  time0 = time2;
728  time1 = sp->time_first;
729  time2 = sp->time_last;
730 
731  // For multi-speed, break out once speed step is beyond data
732  if ( step > 0 && ed->scan[ lscan ].time < sp->time_first )
733  break;
734 
735  // We need to simulate acceleration
736  if ( sp->acceleration_flag )
737  {
738  // If the speed difference is larger than acceleration
739  // rate then we have at least 1 acceleration step
740 
741  af_params.time_steps = (int)
742  ( qAbs( step_speed - current_speed ) / sp->acceleration );
743 
744  // Each simulation step is 1 second long in the acceleration phase
745  af_params.dt = 1.0;
746  af_params.simpoints = 2 * simparams.simpoints;
747 
748  // Use a fixed grid with refinement at both ends and with twice
749  // the number of points
750  duration = time2 - time0;
751  current_time = time0;
752  af_params.start_time = current_time;
753 
754  calculate_ra2( current_speed, step_speed,
755  vC0, simdata, true );
756 int nr5=ed->radius.size();
757 DbgLv(2) << "RSA:(5) nr5" << nr4 << "r sme"
758  << ed->radius[0] << ed->radius[1] << ed->radius[2]
759  << ed->radius[nr5/2-1] << ed->radius[nr5/2] << ed->radius[nr5/2+1]
760  << ed->radius[nr5-3] << ed->radius[nr5-2] << ed->radius[nr5-1];
761 
762  // Add the acceleration time:
763  accel_time = af_params.dt * af_params.time_steps;
764  current_time += accel_time;
765 
766 #ifndef NO_DB
767  emit new_time( current_time );
768  qApp->processEvents();
769 #endif
770 
771  qApp->processEvents();
772  if ( stopFlag ) return 1;
773  } // End of for acceleration
774 
775  duration = time2 - time0;
776  delay = time1 - time0;
777  //duration = ( sp->duration_hours * 3600.
778  // + sp->duration_minutes * 60. );
779 
780  if ( step == simparams.speed_step.size() - 1 )
781  duration += (double)( (int)( duration * 0.05 ) ); // +5%
782 
783  if ( accel_time > duration )
784  {
785  DbgErr() << "Attention: acceleration time exceeds duration - "
786  "please check initialization\n";
787  return -1;
788  }
789  else
790  {
791  duration -= accel_time;
792  }
793  double s_max = qAbs( af_params.s[ 0 ] ); // Find the largest s
794 
795  for ( int mm = 1; mm < af_params.s.size(); mm++ )
796  s_max = qMax( s_max, qAbs( af_params.s[ mm ] ) );
797 
798  af_params.omega_s = sq( step_speed * M_PI / 30. );
799 DbgLv(2) << "RSA: ***COMPUTED(2) omega_s:" << af_params.omega_s << "step_speed"
800  << step_speed << "s.size s_max" << af_params.s.size() << s_max;
801 
802  double lg_bm_rat = log( af_params.current_bottom
803  / af_params.current_meniscus );
804  double s_omg_fac = af_params.omega_s * s_max;
805  af_params.dt = lg_bm_rat
806  / ( s_omg_fac * ( simparams.simpoints - 1 ) );
807 
808  if ( af_params.dt > duration )
809  {
810 double dtsv=af_params.dt;
811  af_params.dt = duration;
812  af_params.simpoints = 1 +
813  (int)( lg_bm_rat / ( s_omg_fac * af_params.dt ) );
814 DbgLv(1) << "RSA: ***CORRECTED(2) dt:" << duration << dtsv << af_params.simpoints;
815  }
816 
817  if ( af_params.simpoints > 10000 )
818  {
819 DbgLv(1) << "RSA: ***COMPUTED simpoints:" << af_params.simpoints
820  << "omg2t s_max dt dura" << af_params.omega_s << s_max << af_params.dt
821  << duration << "** simpts set to 10000 **";
822  af_params.simpoints = 10000;
823  }
824 else
825 DbgLv(2) << "RSA: ***COMPUTED simpoints:" << af_params.simpoints;
826 
827  // Find out the minimum number of simpoints needed to provide the
828  // necessary dt:
829  af_params.time_steps = (int)( duration / af_params.dt ) + 1;
830  af_params.start_time = current_time;
831 US_AstfemMath::MfemData* ed = &af_data;
832 int nr1=ed->radius.size();
833 DbgLv(2) << "RSA:(1) nr1" << nr1 << "r sme"
834  << ed->radius[0] << ed->radius[1] << ed->radius[2]
835  << ed->radius[nr1/2-1] << ed->radius[nr1/2] << ed->radius[nr1/2+1]
836  << ed->radius[nr1-3] << ed->radius[nr1-2] << ed->radius[nr1-1];
837 
838 DbgLv(2) << "RSA: tsteps sttime" << af_params.time_steps << current_time;
839  calculate_ra2( step_speed, step_speed, vC0, simdata, false );
840 
841  // Set the current time to the last scan of this speed step
842  duration = sp->duration_hours * 3600.0
843  + sp->duration_minutes * 60.0;
844  delay = sp->delay_hours * 3600.0
845  + sp->delay_minutes * 60.0;
846  current_time = time1;
847  //current_time = duration;
848  //current_time = time2;
849 DbgLv(2) << "RSA: current_time" << current_time << "fscan lscan"
850  << fscan << lscan << "speed" << step_speed;
851 
852  // Interpolate the simulated data onto the experimental
853  // time and radius grid
854 // if ( ed->scan[ fscan ].time <= sp->time_last )
855 int nr2=ed->radius.size();
856 DbgLv(2) << "RSA:(2) nr2" << nr2 << "r sme"
857  << ed->radius[0] << ed->radius[1] << ed->radius[2]
858  << ed->radius[nr2/2-1] << ed->radius[nr2/2] << ed->radius[nr2/2+1]
859  << ed->radius[nr2-3] << ed->radius[nr2-2] << ed->radius[nr2-1];
860  if ( ed->scan[ fscan ].time <= sp->time_last &&
861  ed->scan[ lscan ].time >= sp->time_first )
862  {
863 // US_AstfemMath::interpolate( *ed, simdata, use_time,
864 // fscan, ++lscan );
865  US_AstfemMath::interpolate( *ed, simdata, use_time );
866  }
867 int nr3=ed->radius.size();
868 DbgLv(2) << "RSA:(3) nr3" << nr3 << "r sme"
869  << ed->radius[0] << ed->radius[1] << ed->radius[2]
870  << ed->radius[nr3/2-1] << ed->radius[nr3/2] << ed->radius[nr3/2+1]
871  << ed->radius[nr3-3] << ed->radius[nr3-2] << ed->radius[nr3-1];
872 
873  if ( !simout_flag )
874  {
875  simdata.radius.clear();
876  simdata.scan .clear();
877  }
878 
879  // Set the current speed to the constant rotor speed of the
880  // current speed step
881  current_speed = step_speed;
882 
883 #ifndef NO_DB
884  qApp->processEvents();
885 #endif
886  if ( stopFlag ) return 1;
887  } // Speed step loop
888 #ifdef TIMING_RA
889 totT5+=(clcSt5.msecsTo(QDateTime::currentDateTime()));
890 #endif
891  } // RG Group loop
892 DbgLv(2) << "RSA: Speed step OUT";
893 #ifdef TIMING_RA
894 QDateTime clcSt7 = QDateTime::currentDateTime();
895 totT6+=(clcSt6.msecsTo(clcSt7));
896 #endif
897 
898 #ifndef NO_DB
899  //emit current_component( -1 );
900  qApp->processEvents();
901 #endif
902 
903 US_AstfemMath::MfemData* ed = &af_data;
904 int nr7=ed->radius.size();
905 DbgLv(2) << "RSA:(7) nr7" << nr7 << "r sme"
906  << ed->radius[0] << ed->radius[1] << ed->radius[2]
907  << ed->radius[nr7/2-1] << ed->radius[nr7/2] << ed->radius[nr7/2+1]
908  << ed->radius[nr7-3] << ed->radius[nr7-2] << ed->radius[nr7-1];
909  if ( time_correction && !simout_flag )
910  {
911  int nstep = simparams.speed_step.size();
912  double correction = 0.0;
913 
914  // Check each speed step to see if it contains acceleration
915  for ( int step = 0; step < nstep; step++ )
916  {
917  US_SimulationParameters::SpeedProfile*
918  sp = &simparams.speed_step[ step ];
919  US_AstfemMath::MfemData* ed = &af_data;
920  int nscans = ed->scan.size();
921 DbgLv(2) << "RSA: TimeCorr step" << step << "nscans" << nscans;
922 
923  if ( nstep > 1 )
924  { // For multi-speed, skip steps outside experimental scan range
925  int fscan = 0;
926  int lscan = nscans - 1;
927 DbgLv(2) << "RSA: TimeCorr nscans fscan lscan" << nscans << fscan << lscan;
928 DbgLv(2) << "RSA: TimeCorr stimef stimel" << sp->time_first << sp->time_last
929  << "etimef etimel" << ed->scan[ fscan ].time << ed->scan[ lscan ].time;
930 
931  if ( ed->scan[ fscan ].time > sp->time_last )
932  continue;
933 
934  if ( step > 0 && ed->scan[ lscan ].time < sp->time_first )
935  break;
936  }
937 
938  // We need to correct time
939  if ( sp->acceleration_flag )
940  {
941  double slope;
942  double intercept;
943  double correlation;
944  double sigma;
945 
946  if ( nscans > 1 )
947  {
948  QVector< double > xtmpVec( nscans );
949  QVector< double > ytmpVec( nscans );
950  double* xtmp = xtmpVec.data();
951  double* ytmp = ytmpVec.data();
952 
953  // Only fit the scans that belong to this speed step
954  for ( int ii = 0; ii < nscans; ii++ )
955  {
956  xtmp[ ii ] = ed->scan[ ii ].time;
957  ytmp[ ii ] = ed->scan[ ii ].omega_s_t;
958  }
959 
960  US_Math2::linefit( &xtmp, &ytmp, &slope, &intercept, &sigma,
961  &correlation, nscans );
962 
963  correction = -intercept / slope;
964  }
965 DbgLv(2) << "RSA: TimeCorr correction" << correction << "nscans" << nscans;
966 
967  for ( int ii = 0; ii < nscans; ii++ )
968  ed->scan[ ii ].time -= correction;
969 
970  // We only need make correction for one speed step
971  break;
972  } // END: ( sp->acceleration_flag )
973 
974  if ( correction > 0.0 )
975  break;
976  } // END: step loop
977  } // END: ( time_correction && !simout_flag )
978 DbgLv(2) << "RSA: Time Corr OUT";
979 #ifdef TIMING_RA
980 QDateTime clcSt8 = QDateTime::currentDateTime();
981 totT7+=(clcSt7.msecsTo(clcSt8));
982 #endif
983 int nr9=ed->radius.size();
984 DbgLv(2) << "RSA:(9) nr9" << nr9 << "r sme"
985  << ed->radius[0] << ed->radius[1] << ed->radius[2]
986  << ed->radius[nr9/2-1] << ed->radius[nr9/2] << ed->radius[nr9/2+1]
987  << ed->radius[nr9-3] << ed->radius[nr9-2] << ed->radius[nr9-1];
988 
989  if ( !simout_flag )
990  store_mfem_data( exp_data, af_data ); // normal experiment grid
991  else
992  store_mfem_data( exp_data, simdata ); // raw simulation grid
993 
994 #ifdef TIMING_RA
995 QDateTime clcSt9 = QDateTime::currentDateTime();
996 totT8+=(clcSt8.msecsTo(clcSt9));
997 int elapsedCalc = calcStart.msecsTo( clcSt9 );
998 ncalls++;
999 totTC+=elapsedCalc;
1000 if((ncalls%TIMING_RA_INC)<1) {
1001  DbgLv(0) << " Elapsed fem-calc ms" << elapsedCalc << "nc totC" << ncalls << totTC << " size_cv" << size_cv;
1002  DbgLv(0) << " t1 t2 t3 t4 t5 t6 t7 t8"
1003  << totT1 << totT2 << totT3 << totT4 << totT5 << totT6 << totT7 << totT8;
1004 }
1005 #endif
1006 DbgLv(2) << "ASTFEM CALC DONE";
1007  return 0;
1008 }
1009 
1010 
1011 void US_Astfem_RSA::update_assocv( void )
1012 {
1013  for ( int ii = 0; ii < system.associations.size(); ii++ )
1014  {
1015  US_Model::Association* as = &system.associations[ ii ];
1016  int ncomp = as->rcomps.size();
1017  int stoich1 = as->stoichs[ 0 ];
1018  int stoich2 = as->stoichs[ 1 ];
1019 DbgLv(2) << "AFRSA: ncomp st1 st2" << ncomp << stoich1 << stoich2;
1020 
1021  if ( ncomp == 2 )
1022  {
1023  as->stoichs[ 1 ] = stoich2 < 0 ? stoich2 : -stoich2;
1024  as->stoichs[ 0 ] = stoich1 > 0 ? stoich1 : -stoich1;
1025  //as->rcomps[ 0 ] = 0;
1026  //as->rcomps[ 1 ] = 1;
1027  }
1028 
1029  else if ( ncomp == 3 )
1030  {
1031  int stoich3 = as->stoichs[ 2 ];
1032  as->stoichs[ 0 ] = stoich1 > 0 ? stoich1 : -stoich1;
1033  as->stoichs[ 1 ] = stoich2 > 0 ? stoich2 : -stoich2;
1034  as->stoichs[ 2 ] = stoich3 < 0 ? stoich3 : -stoich3;
1035  //as->rcomps[ 0 ] = 1;
1036  //as->rcomps[ 1 ] = 1;
1037  //as->rcomps[ 2 ] = 0;
1038  }
1039  }
1040 }
1041 
1042 // Adjust meniscus and bottom based on rotor coefficients
1043 void US_Astfem_RSA::adjust_limits( int speed )
1044 {
1045  // First correct meniscus to theoretical position at rest:
1046  double stretch_value = stretch( simparams.rotorcoeffs,
1047  af_params.first_speed );
1048 
1049  // This is the meniscus at rest
1050  af_params.current_meniscus = simparams.meniscus - stretch_value;
1051 
1052  // Calculate rotor stretch at current speed
1053  stretch_value = stretch( simparams.rotorcoeffs, speed );
1054 
1055  // Add current stretch to meniscus at rest
1056  af_params.current_meniscus += stretch_value;
1057 
1058  // Add current stretch to bottom at rest
1059 // af_params.current_bottom = simparams.bottom + stretch_value;
1060  af_params.current_bottom = simparams.bottom_position + stretch_value;
1061 }
1062 
1063 // Calculate stretch for rotor coefficients array and rpm
1064 double US_Astfem_RSA::stretch( double* rotorcoeffs, int rpm )
1065 {
1066  double speed = (double)rpm;
1067  return ( rotorcoeffs[ 0 ] * speed
1068  + rotorcoeffs[ 1 ] * sq( speed ) );
1069 }
1070 
1071 // Setup reaction groups
1072 void US_Astfem_RSA::initialize_rg( void )
1073 {
1074  rg.clear();
1075 
1076  // If there are no reactions, then it is all noninteracting
1077  if ( system.associations.size() == 0 ) return;
1078  QVector< bool > reaction_used;
1079  reaction_used.clear();
1080 
1081  for ( int i = 0; i < system.associations.size(); i++ )
1082  reaction_used.append( false );
1083 
1084  // Initialize the first reaction group and put it into a temporary reaction
1085  // group, use as test against all assoc vector entries:
1086 
1087  US_Model::Association* as = &system.associations[ 0 ];
1088  US_AstfemMath::ReactionGroup tmp_rg;
1089 
1090  tmp_rg.GroupComponent.clear();
1091  tmp_rg.association.clear();
1092 
1093  tmp_rg.association .append( 0 );
1094  tmp_rg.GroupComponent.append( as->rcomps[ 0 ] );
1095  tmp_rg.GroupComponent.append( as->rcomps[ 1 ] );
1096 
1097  // Only 2 components react in first reaction
1098  if ( as->rcomps.size() > 2 )
1099  tmp_rg.GroupComponent .append( as->rcomps[ 2 ] );
1100 
1101  reaction_used[ 0 ] = true;
1102 
1103  // There is only one reaction, so add it and return
1104  if ( system.associations.size() == 1 )
1105  {
1106  rg .append( tmp_rg );
1107  return;
1108  }
1109 
1110  bool flag3 = false;
1111 
1112  for ( int i = 0; i < system.associations.size(); i++ )
1113  {
1114  int ncomp;
1115  int component1;
1116  int component2;
1117  int component3;
1118 
1119  // Check each association rule to see if it contains components that
1120  // match tmp_rg components
1121 
1122  for ( int counter = 1; counter < system.associations.size(); counter++ )
1123  {
1124  US_Model::Association* av = &system.associations[ counter ];
1125  ncomp = av->rcomps.size();
1126  component1 = ( ncomp > 0 ) ? av->rcomps[ 0 ] : -1;
1127  component2 = ( ncomp > 1 ) ? av->rcomps[ 1 ] : -1;
1128  component3 = ( ncomp > 2 ) ? av->rcomps[ 2 ] : -1;
1129 
1130  while ( reaction_used[ counter ] )
1131  {
1132  counter++;
1133  if ( counter == system.associations.size() ) return;
1134  }
1135 
1136  // Check if any component already present in tmp_rg matches any of the
1137  // three components in current (*system).associations entry
1138 
1139  bool flag1 = false;
1140 
1141  for ( int j = 0; j < tmp_rg.GroupComponent.size(); j++ )
1142  {
1143  flag1 = false;
1144 
1145  if ( component1 == (int) tmp_rg.GroupComponent[ j ] ||
1146  component2 == (int) tmp_rg.GroupComponent[ j ] ||
1147  component3 == (int) tmp_rg.GroupComponent[ j ] )
1148  {
1149  flag1 = true;
1150  break;
1151  }
1152  }
1153 
1154  // If the component from tmp_rg is present in another
1155  // system.associations entry, find out if the component from
1156  // system.associations is already in tmp_rg.GroupComponent
1157 
1158  if ( flag1 )
1159  {
1160  // It is present (flag1=true) so add this rule to the tmp_rg vector
1161  tmp_rg.association .append( counter );
1162  reaction_used[ counter ] = true;
1163 
1164  // There is at least one of all system.associations entries in
1165  // this reaction_group, so set flag3 to true
1166 
1167  flag3 = true;
1168  bool flag2 = false;
1169 
1170  // Check if 1st component is already in the GroupVector from tmp_rg
1171  for ( int j = 0; j < tmp_rg.GroupComponent.size(); j++ )
1172  {
1173  if ( component1 == (int) tmp_rg.GroupComponent[ j ])
1174  flag2 = true;
1175  }
1176 
1177  // Add if not present already
1178  if ( ! flag2 ) tmp_rg.GroupComponent.append( component1 );
1179 
1180  flag2 = false;
1181 
1182  // Check if 2nd component is already in the GroupVector from tmp_rg
1183  for ( int j = 0; j < tmp_rg.GroupComponent.size(); j++ )
1184  {
1185  if ( component2 == (int) tmp_rg.GroupComponent[ j ])
1186  flag2 = true;
1187  }
1188 
1189  // Add if not present already
1190  if ( ! flag2 ) tmp_rg.GroupComponent .append( component2 );
1191 
1192  flag2 = false;
1193 
1194  // Check if 3rd component is already in the GroupVector from tmp_rg
1195  // (but only if non-zero)
1196 
1197  if ( ncomp > 2 )
1198  {
1199  for ( int j = 0; j < tmp_rg.GroupComponent.size(); j++ )
1200  {
1201  if ( component3 == (int) tmp_rg.GroupComponent[ j ] )
1202  flag2 = true;
1203  }
1204 
1205  // Add if not present already
1206  if ( ! flag2 ) tmp_rg.GroupComponent .append( component3 );
1207  }
1208  }
1209  }
1210 
1211  if ( flag3 )
1212  {
1213  flag3 = false;
1214  rg .append( tmp_rg );
1215 
1216  tmp_rg.GroupComponent.clear();
1217  tmp_rg.association.clear();
1218 
1219  // Make the next unused reaction the test reaction
1220  int j = 1;
1221 
1222  while ( reaction_used[ j ] )
1223  {
1224  j++;
1225  if ( j >= reaction_used.size() ) return;
1226  }
1227 
1228  US_Model::Association* avj = &system.associations[ j ];
1229  ncomp = avj->rcomps.size();
1230  component1 = ( ncomp > 0 ) ? avj->rcomps[ 0 ] : 0;
1231  component2 = ( ncomp > 1 ) ? avj->rcomps[ 1 ] : 0;
1232  component3 = ( ncomp > 2 ) ? avj->rcomps[ 2 ] : 0;
1233 
1234  if ( j < system.associations.size() )
1235  {
1236  tmp_rg.association .append( j );
1237  tmp_rg.GroupComponent.append( component1 );
1238  tmp_rg.GroupComponent.append( component2 );
1239 
1240  // Only 2 components react in first reaction
1241  if ( ncomp > 2 )
1242  tmp_rg.GroupComponent.append( component3 );
1243 
1244  reaction_used[ j ] = true;
1245  //counter++; // Out of scope!!!!
1246  }
1247 
1248  if ( j == system.associations.size() - 1 )
1249  {
1250  rg .append( tmp_rg );
1251  return;
1252  }
1253  }
1254  }
1255 }
1256 
1257 // Initializes total concentration vector
1258 void US_Astfem_RSA::initialize_conc( int kk, US_AstfemMath::MfemInitial& CT0,
1259  bool noninteracting )
1260 {
1261 DbgLv(2) << "RSA: init_conc() ENTER kk" << kk
1262  << "af_c0.concsz" << af_c0.concentration.size();
1263  US_Model::SimulationComponent* sc = &system.components[ kk ];
1264  int nval = af_c0.concentration.size();
1265 
1266  // We don't have an existing CT0 concentration vector. Build up the initial
1267  // concentration vector with constant concentration
1268 
1269  if ( nval == 0 )
1270  {
1271 double mxct=0.0;
1272 int jmxc=0;
1273  nval = CT0.concentration.size();
1274  if ( simparams.band_forming )
1275  {
1276  // Calculate the width of the lamella
1277  double angl = simparams.cp_angle != 0.0 ? simparams.cp_angle : 2.5;
1278  double plen = simparams.cp_pathlen != 0.0 ? simparams.cp_pathlen : 1.2;
1279 DbgLv(2) << "RSA: angle pathlen" << angl << plen;
1280 DbgLv(2) << "RSA: bandvol" << simparams.band_volume << " CT0concsz" << CT0.concentration.size();
1281  double base = sq( af_params.current_meniscus )
1282  + simparams.band_volume * 360.0 / ( angl * plen * M_PI );
1283 
1284  double lamella_width = sqrt( base ) - af_params.current_meniscus;
1285 DbgLv(2) << "RSA: menisc base lwid" << af_params.current_meniscus << base << lamella_width;
1286 
1287  // Calculate the spread of the lamella:
1288  for ( int j = 0; j < nval; j++ )
1289  {
1290  base = ( CT0.radius[ j ] - af_params.current_meniscus )
1291  / lamella_width;
1292 
1293  CT0.concentration[ j ] +=
1294  sc->signal_concentration * exp( -pow( base, 4.0 ) );
1295 if(j<2||j>(CT0.concentration.size()-3)||j==(CT0.concentration.size()/40))
1296 DbgLv(2) << "RSA: j base conc" << j << base << CT0.concentration[j];
1297 if(mxct<CT0.concentration[j]) {mxct=CT0.concentration[j];jmxc=j;}
1298  }
1299  }
1300 
1301  else // !simparams.band_forming
1302  {
1303  for ( int j = 0; j < nval; j++ )
1304  {
1305  CT0.concentration[j] += sc->signal_concentration;
1306 if(mxct<CT0.concentration[j]) {mxct=CT0.concentration[j];jmxc=j;}
1307  }
1308  }
1309 DbgLv(2) << "RSA: kk jmxc" << kk << jmxc << "max_conc" << mxct
1310  << "CT0concsz" << CT0.concentration.size();
1311  }
1312 
1313  else // af_c0.concentration.size() > 0
1314  {
1315  if ( noninteracting )
1316  {
1317  // Take the existing initial concentration vector and copy it to the
1318  // temporary CT0 vector: needs rubber band to make sure meniscus and
1319  // bottom equal current_meniscus and current_bottom
1320 
1321  CT0.radius .clear();
1322  CT0.concentration.clear();
1323  CT0.radius .reserve( nval );
1324  CT0.concentration.reserve( nval );
1325  //CT0 = system.components[ k ].c0;
1326  double conc0 = 0.0;
1327  for ( int jj = 0; jj < nval; jj++ )
1328  {
1329  CT0.radius .append( af_c0.radius [ jj ] );
1330  CT0.concentration.append( af_c0.concentration[ jj ] );
1331  conc0 += af_c0.concentration[ jj ];
1332  }
1333  }
1334 
1335  else // interpolation
1336  {
1337  US_AstfemMath::MfemInitial C;
1338  int nval = CT0.concentration.size();
1339  C.radius .clear();
1340  C.concentration.clear();
1341  C.radius .reserve( nval );
1342  C.concentration.reserve( nval );
1343 
1344  double dr = ( af_params.current_bottom - af_params.current_meniscus )
1345  / (double)( nval - 1 );
1346  double rad = af_params.current_meniscus;
1347 
1348  for ( int j = 0; j < nval; j++ )
1349  {
1350  C.radius .append( rad );
1351  C.concentration.append( 0.0 );
1352  rad += dr;
1353  }
1354 
1355  US_AstfemMath::interpolate_C0( af_c0, C );
1356 
1357  for ( int j = 0; j < nval; j++ )
1358  CT0.concentration[ j ] += C.concentration[ j ];
1359  }
1360  }
1361 DbgLv(2) << "RSA: init_conc() RETURN CT0[0] CT[n]" << CT0.concentration[0]
1362  << CT0.concentration[nval-1] << "nval" << nval;
1363 }
1364 
1365 // Non-interacting solute, constant speed
1366 int US_Astfem_RSA::calculate_ni( double rpm_start, double rpm_stop,
1367  US_AstfemMath::MfemInitial& C_init, US_AstfemMath::MfemData& simdata,
1368  bool accel )
1369 {
1370 #ifdef NO_DB
1371  static int Nsave = 0;
1372  static int Nsavea = 0;
1373  static double** CA = NULL; // stiffness matrix on left hand side
1374  // CA[0...Ms-1][0...N-1][4]
1375 
1376  static double** CB = NULL; // stiffness matrix on right hand side
1377  // CB[0...Ms-1][0...N-1][4]
1378 
1379  static double** CA1; // for matrices used in acceleration
1380  static double** CA2;
1381  static double** CB1;
1382  static double** CB2;
1383 #else
1384  double** CA = NULL; // stiffness matrix on left hand side
1385  // CA[0...Ms-1][0...N-1][4]
1386 
1387  double** CB = NULL; // stiffness matrix on right hand side
1388  // CB[0...Ms-1][0...N-1][4]
1389 
1390 
1391  double** CA1; // for matrices used in acceleration
1392  double** CA2;
1393  double** CB1;
1394  double** CB2;
1395 #endif
1396 
1397  double* C0 = NULL; // C[m][j]: current/next concentration of
1398  // m-th component at x_j
1399  double* C1 = NULL; // C[0...Ms-1][0....N-1]:
1400 
1401 #if 0
1402 #define TIMING_NI 1
1403 #endif
1404 #ifdef TIMING_NI
1405 static int nccall=0;
1406 static int ttTC=0;
1407 static int ttT1=0;
1408 static int ttT2=0;
1409 static int ttT3=0;
1410 static int ttT4=0;
1411 static int ttT5=0;
1412 static int ttT6=0;
1413 static int ttT7=0;
1414 static int ttT8=0;
1415 QDateTime clcSt0 = QDateTime::currentDateTime();
1416 QDateTime clcSt1 = clcSt0;
1417 QDateTime clcSt2 = clcSt0;
1418 QDateTime clcSt3 = clcSt0;
1419 QDateTime clcSt4 = clcSt0;
1420 QDateTime clcSt5 = clcSt0;
1421 QDateTime clcSt6 = clcSt0;
1422 QDateTime clcSt7 = clcSt0;
1423 QDateTime clcSt8 = clcSt0;
1424 QDateTime clcSt9 = clcSt0;
1425 #endif
1426 
1427  simdata.radius.clear();
1428  simdata.scan .clear();
1429 DbgLv(2) << "RSA: (4)AP.start_om2t" << af_params.start_om2t;
1430  w2t_integral = af_params.start_om2t;
1431  last_time = af_params.start_time;
1432 
1433  US_AstfemMath::MfemScan simscan;
1434 
1435  // Generate the adaptive mesh
1436 
1437  xA = x.data();
1438  double sw2 = af_params.s[ 0 ] * sq( rpm_stop * M_PI / 30 );
1439  QVector< double > nu;
1440  nu.clear();
1441  nu .append( sw2 / af_params.D[ 0 ] );
1442 
1443  mesh_gen( nu, simparams.meshType );
1444 DbgLv(2) << "RSA: mesh_gen size" << nu.count() << "accel" << accel;
1445 
1446  // Refine left hand side (when s > 0) or right hand side (when s < 0)
1447  // for acceleration
1448 
1449  if ( accel )
1450  {
1451  int j;
1452  double xc;
1453 
1454  if ( af_params.s[ 0 ] > 0 )
1455  {
1456  // Radial distance from meniscus how far the boundary will move during
1457  // this acceleration step (without diffusion)
1458 
1459  xc = af_params.current_meniscus +
1460  sw2 * ( af_params.time_steps * af_params.dt ) / 3.0;
1461 
1462  for ( j = 0; j < Nx - 3; j++ )
1463  if ( xA[ j ] > xc ) break;
1464  }
1465  else
1466  {
1467  xc = af_params.current_bottom +
1468  sw2 * ( af_params.time_steps * af_params.dt ) / 3.0;
1469 
1470  for ( j = 0; j < Nx - 3; j++ )
1471  if ( xA[ Nx - j - 1 ] < xc ) break;
1472  }
1473 
1474  mesh_gen_RefL( j + 1, 4 * j );
1475  }
1476 
1477 //DbgLv(2) << "RSA: simdat radsize scnsize" << simdata.radius.size()
1478 // << simdata.scan.size() << "Nx" << Nx;
1479  simdata.radius.clear();
1480  simdata.scan .clear();
1481  simdata.radius.reserve( Nx );
1482 // simdata.scan .reserve( Nx );
1483 
1484  for ( int i = 0; i < Nx; i++ )
1485  simdata.radius .append( xA[ i ] );
1486 
1487  // Initialize the coefficient matrices
1488 
1489 #ifdef NO_DB
1490  if ( Nx > Nsave )
1491  {
1492  if ( Nsave > 0 )
1493  {
1494  US_AstfemMath::clear_2d( 3, CA );
1495  US_AstfemMath::clear_2d( 3, CB );
1496  }
1497 
1498  Nsave = Nx;
1499  US_AstfemMath::initialize_2d( 3, Nsave, &CA );
1500  US_AstfemMath::initialize_2d( 3, Nsave, &CB );
1501  }
1502  else
1503  {
1504  for ( int ii = 0; ii < 3; ii++ )
1505  {
1506  for ( int jj = 0; jj < Nx; jj++ )
1507  {
1508  CA[ ii ][ jj ] = 0.0;
1509  CB[ ii ][ jj ] = 0.0;
1510  }
1511  }
1512  }
1513 #else
1514  US_AstfemMath::initialize_2d( 3, Nx, &CA );
1515  US_AstfemMath::initialize_2d( 3, Nx, &CB );
1516 #endif
1517 
1518  bool fixedGrid = ( simparams.gridType == US_SimulationParameters::FIXED );
1519 #ifdef TIMING_NI
1520 clcSt2 = QDateTime::currentDateTime();
1521 ttT1+=(clcSt1.msecsTo(clcSt2));
1522 #endif
1523 
1524  if ( ! accel ) // No acceleration
1525  {
1526  sw2 = af_params.s[ 0 ] * sq( rpm_stop * M_PI / 30 );
1527 
1528 //DbgLv(2) << "RSA: fixedGrid s0" << fixedGrid << af_params.s[0];
1529  if ( fixedGrid )
1530  {
1531  ComputeCoefMatrixFixedMesh ( af_params.D[ 0 ], sw2, CA, CB );
1532  }
1533  else if ( af_params.s[ 0 ] > 0 )
1534  {
1535  ComputeCoefMatrixMovingMeshR( af_params.D[ 0 ], sw2, CA, CB );
1536  }
1537  else
1538  {
1539  ComputeCoefMatrixMovingMeshL( af_params.D[ 0 ], sw2, CA, CB );
1540  }
1541  }
1542  else // For acceleration
1543  {
1544 #ifdef NO_DB
1545  if ( Nx > Nsavea )
1546  {
1547  if ( Nsavea > 0 )
1548  {
1549  US_AstfemMath::clear_2d( 3, CA1 );
1550  US_AstfemMath::clear_2d( 3, CB1 );
1551  US_AstfemMath::clear_2d( 3, CA2 );
1552  US_AstfemMath::clear_2d( 3, CB2 );
1553  }
1554 
1555  US_AstfemMath::initialize_2d( 3, Nx, &CA1 );
1556  US_AstfemMath::initialize_2d( 3, Nx, &CA2 );
1557  US_AstfemMath::initialize_2d( 3, Nx, &CB1 );
1558  US_AstfemMath::initialize_2d( 3, Nx, &CB2 );
1559  Nsavea = Nx;
1560  }
1561  else
1562  {
1563  for ( int ii = 0; ii < 3; ii++ )
1564  {
1565  for ( int jj = 0; jj < Nx; jj++ )
1566  {
1567  CA1[ ii ][ jj ] = 0.0;
1568  CA2[ ii ][ jj ] = 0.0;
1569  CB1[ ii ][ jj ] = 0.0;
1570  CB2[ ii ][ jj ] = 0.0;
1571  }
1572  }
1573  }
1574 #else
1575  US_AstfemMath::initialize_2d( 3, Nx, &CA1 );
1576  US_AstfemMath::initialize_2d( 3, Nx, &CA2 );
1577  US_AstfemMath::initialize_2d( 3, Nx, &CB1 );
1578  US_AstfemMath::initialize_2d( 3, Nx, &CB2 );
1579 #endif
1580 
1581  sw2 = 0.0;
1582  ComputeCoefMatrixFixedMesh( af_params.D[ 0 ], sw2, CA1, CB1 );
1583 
1584  sw2 = af_params.s[ 0 ] * sq( rpm_stop * M_PI / 30 );
1585  ComputeCoefMatrixFixedMesh( af_params.D[ 0 ], sw2, CA2, CB2 );
1586  }
1587 #ifdef TIMING_NI
1588 clcSt3 = QDateTime::currentDateTime();
1589 ttT2+=(clcSt2.msecsTo(clcSt3));
1590 #endif
1591 
1592  // Initial condition
1593  QVector< double > C0Vec( Nx );
1594  QVector< double > C1Vec( Nx );
1595  QVector< double > rhVec( Nx );
1596  C0 = C0Vec.data();
1597  C1 = C1Vec.data();
1598 
1599  // Interpolate the given C_init vector on the new C0 grid
1600 //DbgLv(2) << "RSA: interp C0";
1601  US_AstfemMath::interpolate_C0( C_init, C0, x );
1602 //DbgLv(2) << "RSA: interp C0 RTN";
1603 
1604  // Time evolution
1605  double* right_hand_side = rhVec.data();
1606 #ifdef TIMING_NI
1607 ttT3+=(clcSt3.msecsTo(QDateTime::currentDateTime()));
1608 clcSt3 = QDateTime::currentDateTime();
1609 #endif
1610  double time_dif = ( af_params.time_steps > 1 && rpm_stop != rpm_start )
1611  ? (double)( af_params.time_steps ) : 1.0;
1612 // double time_dif = ( af_params.time_steps > 2 && rpm_stop != rpm_start )
1613 // ? (double)( af_params.time_steps - 1 ) : 1.0;
1614  double rpm_inc = ( rpm_stop - rpm_start ) / time_dif;
1615  double rpm_current = rpm_start - rpm_inc;
1616  int ntsteps = af_params.time_steps;
1617 // int nisteps = ntsteps + 2;
1618  int nisteps = ntsteps + 3;
1619 // int ltsteps = ntsteps - 1;
1620  int ltsteps = ntsteps;
1621 //DbgLv(2) << "RSA: nisteps" << nisteps;
1622  simdata.scan .reserve( nisteps );
1623 
1624  // Calculate all time steps (plus a little overlap)
1625  for ( int ii = 0; ii < nisteps; ii++ )
1626  {
1627 //DbgLv(2) << "RSA: ii nisteps" << ii << nisteps;
1628 #ifdef TIMING_NI
1629 clcSt4 = QDateTime::currentDateTime();
1630 ttT3+=(clcSt3.msecsTo(clcSt4));
1631 #endif
1632 // if ( ii < ntsteps )
1633  if ( ii <= ntsteps )
1634  rpm_current += rpm_inc;
1635 
1636 #ifndef NO_DB
1637  emit current_speed( (int) rpm_current );
1638 #endif
1639 
1640 //DbgLv(2) << "RSA: (2) ii" << ii << "accel" << accel;
1641  if ( accel ) // Then we have acceleration
1642  {
1643  double rpm_ratio = sq( rpm_current / rpm_stop );
1644 
1645  for ( int j1 = 0; j1 < 3; j1++ )
1646  {
1647  double* CAj = CA [ j1 ];
1648  double* CBj = CB [ j1 ];
1649  double* CA1j = CA1[ j1 ];
1650  double* CA2j = CA2[ j1 ];
1651  double* CB1j = CB1[ j1 ];
1652  double* CB2j = CB2[ j1 ];
1653 
1654  for ( int j2 = 0; j2 < Nx; j2++ )
1655  {
1656  CAj[ j2 ] = CA1j[ j2 ] + rpm_ratio * ( CA2j[ j2 ] - CA1j[ j2 ] );
1657  CBj[ j2 ] = CB1j[ j2 ] + rpm_ratio * ( CB2j[ j2 ] - CB1j[ j2 ] );
1658  }
1659  }
1660  }
1661 //DbgLv(2) << "RSA: (3) ii" << ii;
1662 
1663  simscan.rpm = (int) rpm_current;
1664  simscan.time = af_params.start_time + ii * af_params.dt;
1665 
1666  w2t_integral += ( simscan.time - last_time ) *
1667  sq( rpm_current * M_PI / 30.0 );
1668 if(ii<3||ii>ntsteps-2) {
1669 DbgLv(2) << "TMS:RSA:ni: iistep" << ii << "time ltime omg"
1670  << simscan.time << last_time << w2t_integral << "rpm" << rpm_current
1671  << "st_omg" << af_params.start_om2t << "Nx" << Nx;}
1672 
1673  last_time = simscan.time;
1674  simscan.omega_s_t = w2t_integral;
1675  simscan.temperature = af_data.scan[ 0 ].temperature;
1676 DbgLv(2) << "TMS:RSA:ni: time omega_s_t" << simscan.time << simscan.omega_s_t
1677  << "rpm_c" << rpm_current << "step-scan" << simdata.scan.size();
1678 
1679  simscan.conc.clear();
1680  simscan.conc.reserve( Nx );
1681 DbgLv(2) << "RSA: Nx C0size" << Nx << C0Vec.size() << "step-scan"
1682  << simdata.scan.size() << "rss-now" << US_Memory::rss_now();
1683 
1684  for ( int jj = 0; jj < Nx; jj++ )
1685  simscan.conc.append( C0[ jj ] );
1686 
1687 //DbgLv(2) << "RSA: (4) ii" << ii << "simdata scans"
1688 // << simdata.scan.size();
1689  simdata.scan.append( simscan );
1690 //DbgLv(2) << "RSA: (5) ii" << ii;
1691 if(ii==0) DbgLv(2) << "TMS:RSA:ni: Scan Added";
1692 #ifdef TIMING_NI
1693 clcSt5 = QDateTime::currentDateTime();
1694 ttT4+=(clcSt4.msecsTo(clcSt5));
1695 #endif
1696 
1697  // Sedimentation part:
1698  // Calculate the right hand side vector
1699 
1700  if ( accel || fixedGrid )
1701  {
1702  right_hand_side[ 0 ] = - CB[ 1 ][ 0 ] * C0[ 0 ]
1703  - CB[ 2 ][ 0 ] * C0[ 1 ];
1704 
1705  for ( int jj = 1; jj < Nx - 1; jj++ )
1706  {
1707  right_hand_side[ jj ] = - CB[ 0 ][ jj ] * C0[ jj - 1 ]
1708  - CB[ 1 ][ jj ] * C0[ jj ]
1709  - CB[ 2 ][ jj ] * C0[ jj + 1 ];
1710  }
1711 
1712  int jj = Nx - 1;
1713  right_hand_side[ jj ] = - CB[ 0 ][ jj ] * C0[ jj - 1 ]
1714  - CB[ 1 ][ jj ] * C0[ jj ];
1715  }
1716  else
1717  {
1718  if ( af_params.s[ 0 ] > 0 )
1719  {
1720  right_hand_side[ 0 ] = - CB[ 2 ][ 0 ] * C0[ 0 ];
1721  right_hand_side[ 1 ] = - CB[ 1 ][ 1 ] * C0[ 0 ]
1722  - CB[ 2 ][ 1 ] * C0[ 1 ];
1723 
1724  for ( int jj = 2; jj < Nx; jj++ )
1725  {
1726  right_hand_side[ jj ] = - CB[ 0 ][ jj ] * C0[ jj - 2 ]
1727  - CB[ 1 ][ jj ] * C0[ jj - 1 ]
1728  - CB[ 2 ][ jj ] * C0[ jj ];
1729  }
1730  }
1731  else
1732  {
1733  for ( int jj = 0; jj < Nx - 2; jj++ )
1734  {
1735  right_hand_side[ jj ] = - CB[ 0 ][ jj ] * C0[ jj ]
1736  - CB[ 1 ][ jj ] * C0[ jj + 1 ]
1737  - CB[ 2 ][ jj ] * C0[ jj + 2 ];
1738  }
1739 
1740  int jj = Nx - 2;
1741  right_hand_side[ jj ] = - CB[ 0 ][ jj ] * C0[ jj ]
1742  - CB[ 1 ][ jj ] * C0[ jj + 1 ];
1743 
1744  jj = Nx - 1;
1745  right_hand_side[ jj ] = -CB[ 0 ][ jj ] * C0[ jj ];
1746  }
1747  }
1748 //DbgLv(2) << "RSA: (6) ii" << ii;
1749 
1750 #ifdef TIMING_NI
1751 clcSt6 = QDateTime::currentDateTime();
1752 ttT5+=(clcSt5.msecsTo(clcSt6));
1753 #endif
1754 //DbgLv(2) << "RSA: tridiag";
1755  US_AstfemMath::tridiag( CA[0], CA[1], CA[2], right_hand_side, C1, Nx );
1756 //DbgLv(2) << "RSA: tridiag RTN";
1757 #ifdef TIMING_NI
1758 clcSt7 = QDateTime::currentDateTime();
1759 ttT6+=(clcSt6.msecsTo(clcSt7));
1760 #endif
1761 
1762 //DbgLv(2) << "RSA: C1size C0size" << C1Vec.size() << C0Vec.size();
1763  for ( int jj = 0; jj < Nx; jj++ ) C0[ jj ] = C1[ jj ];
1764 //DbgLv(2) << "RSA: (7) ii" << ii;
1765 
1766 #ifndef NO_DB
1767  if ( show_movie )
1768  {
1769  qApp->processEvents();
1770  if ( stopFlag ) break;
1771 
1772  emit new_scan( &x, C0 );
1773  emit new_time( simscan.time );
1774  qApp->processEvents();
1775  }
1776 #endif
1777 #ifdef TIMING_NI
1778 clcSt3 = QDateTime::currentDateTime();
1779 ttT7+=(clcSt7.msecsTo(clcSt3));
1780 #endif
1781 
1782 //DbgLv(2) << "RSA: ii ltsteps" << ii << ltsteps;
1783  if ( ii == ltsteps )
1784  {
1785  C_init.radius .clear();
1786  C_init.concentration.clear();
1787  C_init.radius .reserve( Nx );
1788  C_init.concentration.reserve( Nx );
1789  for ( int jj = 0; jj < Nx; jj++ )
1790  {
1791  C_init.radius .append( x [ jj ] );
1792  C_init.concentration .append( C1[ jj ] );
1793  }
1794  }
1795  } // time loop
1796 
1797 #ifdef TIMING_NI
1798 clcSt8 = QDateTime::currentDateTime();
1799 #endif
1800 
1801 #ifndef NO_DB
1802  emit new_time( simscan.time );
1803  qApp->processEvents();
1804  US_AstfemMath::clear_2d( 3, CA );
1805  US_AstfemMath::clear_2d( 3, CB );
1806 
1807  if ( accel )
1808  {
1809  US_AstfemMath::clear_2d( 3, CA1 );
1810  US_AstfemMath::clear_2d( 3, CB1 );
1811  US_AstfemMath::clear_2d( 3, CA2 );
1812  US_AstfemMath::clear_2d( 3, CB2 );
1813  }
1814 #endif
1815 //DbgLv(2) << "RSA: calc_ni() RETURN Nx" << Nx << "C_init[0] C_init[n]"
1816 // << C_init.concentration[0] << C_init.concentration[Nx-1];
1817 
1818 #ifdef TIMING_NI
1819 clcSt9 = QDateTime::currentDateTime();
1820 ttT8+=(clcSt8.msecsTo(clcSt9));
1821 int elapsCalc=clcSt0.msecsTo(clcSt9);
1822 nccall++;
1823 ttTC+=elapsCalc;
1824 if((nccall%1000)==0) {
1825  int totTK=ttT1+ttT2+ttT3+ttT4+ttT5+ttT6+ttT7+ttT8;
1826  qDebug() << " ++Elapsed calc-ni ms" << elapsCalc << " nc totC totK" << nccall
1827  << ttTC << totTK << " Nx" << Nx << " acc fixG" << accel << fixedGrid;
1828  qDebug() << " CN: t1 t2 t3 t4 t5 t6 t7 t8" << ttT1 << ttT2 << ttT3 << ttT4
1829  << ttT5 << ttT6 << ttT7 << ttT8;
1830  int p0=totTK/2;
1831  int p1=(ttT1*100+p0)/totTK;
1832  int p2=(ttT2*100+p0)/totTK;
1833  int p3=(ttT3*100+p0)/totTK;
1834  int p4=(ttT4*100+p0)/totTK;
1835  int p5=(ttT5*100+p0)/totTK;
1836  int p6=(ttT6*100+p0)/totTK;
1837  int p7=(ttT7*100+p0)/totTK;
1838  int p8=(ttT8*100+p0)/totTK;
1839  qDebug() << " CN: p1 p2 p3 p4 p5 p6 p7 p8" << p1 << p2 << p3 << p4
1840  << p5 << p6 << p7 << p8;
1841 }
1842 #endif
1843 DbgLv(2) << "RSA: CALC_NI END - simdata scans points rss"
1844  << simdata.scan.size() << simdata.scan[0].conc.size() << US_Memory::rss_now();
1845 
1846  return 0;
1847 }
1848 
1849 void US_Astfem_RSA::mesh_gen( QVector< double >& nu, int MeshOpt )
1850 {
1852 //
1853 // Generate adaptive grids for multi-component Lamm equations
1854 //
1855 //
1856 // Here: Nx: Number of points in the ASTFEM
1857 // m, b: meniscus, bottom
1858 // nuMax, nuMin = max and min of nu=sw^2/D
1859 // MeshType: = 0 ASTFEM grid based on all nu (composite in sharp region)
1860 // = 1 Claverie (uniform), etc,
1861 // = 2 Exponential mesh (Schuck's form., no refinement at bottom)
1862 // = 3 input from data file: "mesh_data.dat"
1863 // = 4 ASTFVM grid (Finite Volume Method)
1864 // = 10, acceleration mesh (left and right refinement)
1866 
1868 // generate the mesh
1870 
1871  double m = af_params.current_meniscus;
1872  double b = af_params.current_bottom;
1873  int Np = af_params.simpoints;
1874 
1875  x.clear();
1876  x.reserve( Np * 2 + 2 );
1877 
1878  switch ( MeshOpt )
1879  {
1881  // Mesh Type 0 (default): adaptive mesh based on all nu
1883 
1884  case (int)US_SimulationParameters::ASTFEM:
1885  // Adaptive Space Time FE Mesh without left hand refinement
1886  qSort( nu ); // put nu in ascending order
1887 
1888  if ( nu[ 0 ] > 0 )
1889  mesh_gen_s_pos( nu );
1890 
1891  else if ( nu[ nu.size() - 1 ] < 0 )
1892  mesh_gen_s_neg( nu );
1893 
1894  else // Some species with s < 0 and some with s > 0
1895  {
1896  double bmval = m;
1897  double deltbm = ( b - m ) / (double)( Np - 1 );
1898 
1899  for ( int i = 0; i < Np; i++ )
1900  {
1901  x .append( bmval );
1902  bmval += deltbm;
1903  }
1904  }
1905  break;
1906 
1907  case (int)US_SimulationParameters::CLAVERIE:
1908  // Claverie mesh without left hand refinement
1909 
1910  for ( int i = 0; i < Np; i++ )
1911  x .append( m + ( b - m ) * i / ( Np - 1 ) );
1912  break;
1913 
1914  case (int)US_SimulationParameters::MOVING_HAT:
1915  // Moving Hat (Peter Schuck's Mesh) w/o left hand side refinement
1916 
1917  x .append( m );
1918 
1919  // Standard Schuck grids
1920  for ( int i = 1; i < Np - 1; i++ )
1921  x .append( m * pow( b / m, ( i - 0.5 ) / ( Np - 1 ) ) );
1922 
1923  x .append( b );
1924  break;
1925 
1926  case (int)US_SimulationParameters::USER:
1927  // User defined mesh generated from data file
1928  {
1929  QFile f( US_Settings::appBaseDir() + "/etc/mesh.dat" );
1930 
1931  if ( f.open( QIODevice::ReadOnly ) )
1932  {
1933  QTextStream ts( &f );
1934 
1935  while ( ! ts.atEnd() ) x .append( ts.readLine().toDouble() );
1936 
1937  f.close();
1938 
1939  if ( qAbs( x[ 0 ] - m ) > 1.0e7 )
1940  {
1941  DbgErr() << "The meniscus from the mesh file does not"
1942  " match the set meniscus - using Claverie Mesh instead\n";
1943  }
1944 
1945  if ( qAbs( x[ x.size() - 1 ] - b ) > 1.0e7 )
1946  {
1947  DbgErr() << "The cell bottom from the mesh file does not"
1948  " match the set cell bottom - using Claverie Mesh"
1949  " instead\n";
1950  }
1951  }
1952  else
1953  {
1954  DbgErr() << "Could not read the mesh file - "
1955  "using Claverie Mesh instead\n";
1956 
1957  for ( int i = 0; i < af_params.simpoints; i++ )
1958  x .append( m + ( b - m ) * i / ( Np - 1 ) );
1959  }
1960  break;
1961  }
1962 
1963  case (int)US_SimulationParameters::ASTFVM:
1964  // Adaptive Space Time Finite Volume Method
1965  qSort( nu ); // put nu in ascending order
1966 
1967  if ( nu[ 0 ] > 0 )
1968  mesh_gen_s_pos( nu );
1969 
1970  else if ( nu[ nu.size() - 1 ] < 0 )
1971  mesh_gen_s_neg( nu );
1972 
1973  else // Some species with s < 0 and some with s > 0
1974  {
1975  for ( int i = 0; i < Np; i++ )
1976  x .append( m + ( b - m ) * i / ( Np - 1 ) );
1977  }
1978  break;
1979 
1980  default:
1981  qDebug() << "undefined mesh option\n";
1982  break;
1983 
1984  }
1985 
1986  Nx = x.size();
1987  xA = x.data();
1988 DbgLv(2) << "RSA: ***COMPUTED Nx" << Nx << "simpts" << af_params.simpoints
1989  << "omg2t" << af_params.omega_s << "s0" << af_params.s[0]
1990  << "dt" << af_params.dt;
1991 }
1992 
1994 //
1995 // Generate exponential mesh and refine cell bottom (for s>0)
1996 //
1998 void US_Astfem_RSA::mesh_gen_s_pos( const QVector< double >& nuvec )
1999 {
2000  double tmp_Hstar;
2001  QVector< double > xc;
2002  QVector< double > Hstar;
2003  QVector< double > y;
2004  QVector< int > Nf;
2005  const double* nu = nuvec.data();
2006 
2007  xc .clear();
2008  Hstar.clear();
2009  Nf .clear();
2010  xc .reserve( af_params.s.size() * 2 );
2011  Hstar.reserve( af_params.s.size() );
2012  Nf .reserve( af_params.s.size() );
2013 
2014  double meniscus = af_params.current_meniscus;
2015  double bottom = af_params.current_bottom;
2016  int Np = af_params.simpoints;
2017 
2018  int IndLayer = 0; // number of layers for grids in steep region
2019  double uth = 1.0 / Np; // threshold of u for steep region
2020  double bmsqdf = ( sq( bottom ) - sq( meniscus ) ) / ( uth * 2.0 );
2021  double Npm1 = (double)( Np - 1 );
2022  double Npm2 = Npm1 - 1.0;
2023  double Npratio = Npm2 / Npm1;
2024  double bmratio = bottom / meniscus;
2025  double mbratio = meniscus / bottom;
2026  double bmrlog = log( bmratio );
2027  double k2log = log( 2.0 );
2028  double bmNpow = pow( bmratio, Npratio );
2029  double bmdiff = bottom - meniscus * bmNpow;
2030  const double PIhalf = M_PI / 2.0;
2031 DbgLv(2) << "RSA: msg_pos: nu0" << nu[0] << "nu1" << nu[1] << "Np" << Np << "uth" << uth;
2032 DbgLv(2) << "RSA: msg_pos: m " << meniscus << " b " << bottom << " bmsqdf" << bmsqdf;
2033 
2034  for ( int i = 0; i < af_params.s.size(); i++ ) // Markers for steep regions
2035  {
2036  double tmp_xc = bottom - ( 1.0 / ( nu[ i ] * bottom ) )
2037  * log( nu[ i ] * bmsqdf );
2038 
2039  // # of pts for i-th layer
2040  int tmp_Nf = (int) ( PIhalf * ( bottom - tmp_xc )
2041  * nu[ i ] * bottom / 2.0 + 0.5 ) + 1;
2042 
2043  // Step required by Pac(i) < 1
2044  tmp_Hstar = ( bottom - tmp_xc ) / tmp_Nf * PIhalf;
2045 
2046 DbgLv(2) << "RSA: msg_pos: i " << i << "txc" << tmp_xc << "tNf" << tmp_Nf << "tHs" << tmp_Hstar;
2047  if ( ( tmp_xc > meniscus ) && ( bmdiff > tmp_Hstar ) )
2048  {
2049  xc .append( tmp_xc );
2050  Nf .append( tmp_Nf );
2051  Hstar .append( tmp_Hstar );
2052  IndLayer++;
2053  }
2054  }
2055 
2056  xc .append( bottom );
2057 
2058  int indp = 0; // Index for a grid point
2059 
2060  if ( IndLayer > 0 )
2061  { // Steep region (potentially)
2062 
2063  for ( int i = 0; i < IndLayer; i++ ) // Consider i-th steep region
2064  {
2065  double xci = xc[ i ];
2066 
2067  if ( i < IndLayer - 1 ) // Linear distribution for step size distrib
2068  {
2069  double xcip = xc[ i + 1 ];
2070  double HL = Hstar[ i ];
2071  double HR = Hstar[ i + 1 ];
2072  int Mp = (int) ( ( xcip - xci ) * 2.0 / ( HL + HR ) );
2073 
2074  if ( Mp > 1 )
2075  {
2076  double beta = ( ( HR - HL ) / 2.0 ) * Mp;
2077  double alpha = xcip - xci - beta;
2078 
2079  for ( int j = 0; j <= Mp - 1; j++ )
2080  {
2081  double xi = (double) j / (double) Mp;
2082  y .append( xci + alpha * xi + beta * sq( xi ) );
2083  indp++;
2084  }
2085  }
2086  }
2087  else // Last layer, use sine distribution for grids
2088  {
2089  double xcip = xc[ qMin( i + 1, xc.size() - 1 ) ];
2090 
2091  for ( int j = 0; j <= Nf[ i ] - 1; j++ )
2092  {
2093  indp++;
2094  y .append( xci + ( bottom - xci ) *
2095  sin( j / ( Nf[ i ] - 1.0 ) * PIhalf ) );
2096 
2097  if ( y[ indp - 1 ] > xcip ) break;
2098  }
2099  }
2100  }
2101  }
2102 
2103  if ( indp < 2 )
2104  { // IndLayer==0 or indp count less than 2
2105  x .append( meniscus );
2106 
2107  // Add one more point to Schuck's grids
2108  for ( int k = 1; k < Np - 1 ; k++ )
2109  { // Schuck's mesh
2110  x .append( meniscus * pow( bmratio, (double) k / Npm1 ) );
2111  }
2112 
2113  x .append( bottom );
2114  }
2115 
2116  else
2117  { // IndLayer>0 and indp greater than 1
2118  int NfTotal = indp;
2119 DbgLv(2) << "RSA: msg_pos: IndL>0 x sz" << x.size();
2120 
2121  // Reverse the order of y
2122  int jj = 0;
2123  int kk = NfTotal - 1;
2124  double* yary = y.data();
2125 
2126  while ( jj < kk )
2127  {
2128  double ysav = yary[ jj ];
2129  yary[ jj++ ] = yary[ kk ];
2130  yary[ kk-- ] = ysav;
2131  }
2132 
2133  // Transition region
2134  // Smallest step size in transit region
2135  double Hf = y[ NfTotal - 2 ] - y[ NfTotal - 1 ];
2136 
2137  // Number of pts in trans region
2138  int Nm = (int)( log( bottom / ( Npm1 * Hf ) * bmrlog ) / k2log ) + 1;
2139 
2140  double xa = y[ NfTotal - 1 ] - Hf * ( pow( 2.0, (double)Nm ) - 1.0 );
2141 
2142  int Js = (int)( Npm1 * log( xa / meniscus ) / bmrlog );
2143 
2144  // xa is modified so that y[ NfTotal - Nm ] matches xa exactly
2145  xa = meniscus * pow( bmratio, ( (double)Js / (double)Npm1 ) );
2146 
2147  double tmp_xc = y[ NfTotal - 1 ];
2148  double HL = xa * ( 1.0 - mbratio );
2149  double HR = y[ NfTotal - 2 ] - y[ NfTotal - 1 ];
2150 
2151  int Mp = (int)( ( tmp_xc - xa ) * 2.0 / ( HL + HR ) ) + 1;
2152 
2153  if ( Mp > 1 )
2154  {
2155  double beta = ( ( HR - HL ) / 2.0 ) * Mp;
2156  double alpha = ( tmp_xc - xa ) - beta;
2157 
2158  for ( int j = Mp - 1; j > 0; j-- )
2159  {
2160  double xi = (double) j / Mp;
2161  y .append( xa + alpha * xi + beta * sq( xi ) );
2162  }
2163  }
2164 
2165  Nm = Mp;
2166 DbgLv(2) << "RSA: msg_pos: IndL>0 Js" << Js << "NfTotal" << NfTotal << "Nm" << Nm;
2167 
2168  // Regular region
2169  x .append( meniscus );
2170  yary = y.data();
2171 
2172  for ( int j = 1; j <= Js; j++ )
2173  x .append( meniscus * pow( bmratio, ( (double)j / (double)Npm1 ) ) );
2174 
2175  for ( int j = NfTotal + Nm - 2; j >=0; j-- )
2176  x .append( yary[ j ] );
2177 DbgLv(2) << "RSA: msg_pos: IndL>0 y sz" << y.size() << "Mp" << Mp << "Nm" << Nm << "x sz" << x.size();
2178  }
2179 
2180  xA = x.data();
2181 int mm=x.size()/2;
2182 int ee=x.size()-1;
2183 DbgLv(2) << "RSA: mgs_pos: xA sme" << x[0] << x[1] << x[2]
2184  << x[mm-1] << x[mm] << x[mm+1] << x[mm+2] << x[ee-2] << x[ee-1] << x[ee];
2185 }
2186 
2188 //
2189 // Generate exponential mesh and refine cell meniscus (for s<0)
2190 //
2192 void US_Astfem_RSA::mesh_gen_s_neg( const QVector< double >& nu )
2193 {
2194 
2195  const double PIhalf = M_PI / 2.0;
2196  const double PIquar = M_PI / 4.0;
2197  const double k2log = log( 2.0 );
2198  int j, Js, Nf, Nm;
2199  double xc, xa, Hstar;
2200  QVector< double > yr, ys, yt;
2201 
2202  x .clear();
2203  yr.clear();
2204  ys.clear();
2205  yt.clear();
2206 
2207  double m = af_params.current_meniscus;
2208  double b = af_params.current_bottom;
2209  int Np = af_params.simpoints;
2210 
2211  double uth = 1.0 / Np; // Threshold of u for steep region
2212  double nu0 = qAbs( nu[ 0 ] );
2213  double bmrlog = log( b / m );
2214  double mbratio = m / b;
2215  double Npm1 = (double)( Np - 1 );
2216 
2217  x .reserve( Np );
2218  yr.reserve( Np );
2219  ys.reserve( Np );
2220  yt.reserve( Np );
2221 
2222  xc = m + 1. / ( nu0 * m ) * log( ( sq( b ) - sq( m ) ) * nu0 / ( 2. * uth ) );
2223 
2224  Nf = 1 + (int)( ( xc - m ) * nu0 * m * PIquar );
2225 
2226  Hstar = ( xc - m ) / Nf * PIhalf;
2227 
2228  Nm = 1 + (int)( log( m / ( Npm1 * Hstar ) * bmrlog ) / k2log );
2229 
2230  xa = xc + ( pow( 2.0, (double) Nm ) - 1.0 ) * Hstar;
2231 
2232  Js = (int) ( Npm1 * log( b / xa ) / bmrlog + 0.5 );
2233 
2234 
2235  // All grid points at exponentials
2236  yr .append( b );
2237 
2238  // Is there a difference between simparams.meniscus and
2239  // af_params.current_meniscus??
2240  for( j = 1; j < Np; j++ ) // Add one more point to Schuck's grids
2241  yr .append( b * pow( simparams.meniscus / b, ( j - 0.5 ) / Npm1 ) );
2242 
2243  yr .append( m );
2244 
2245  if ( b * ( pow( mbratio, ( Np - 3.5 ) / Npm1 )
2246  - pow( mbratio, ( Np - 2.5 ) / Npm1 ) ) < Hstar || Nf <= 2 )
2247  {
2248  double* yrA = yr.data();
2249 
2250  // No need for steep region
2251  for ( j = Np - 1; j >= 0; j-- )
2252  {
2253  x .append( yrA[ j ] );
2254  }
2255 
2256  xA = x.data();
2257 
2258  DbgErr() << "Use exponential grid only!(1/10000 reported): Np Nf b m nu0"
2259  << Np << Nf << b << m << nu0;
2260  }
2261  else
2262  {
2263  // Nf > 2
2264  double xcm = xc - m;
2265  double Nfm1 = (double)( Nf - 1 );
2266  for ( j = 0; j < Nf - 1; j++ )
2267  ys .append( xc - xcm * sin( (double)j / Nfm1 * PIhalf ) );
2268 
2269  ys .append( m );
2270 
2271  for ( j = 0; j < Nm; j++ )
2272  yt .append( xc + ( pow( 2.0, (double) j ) - 1.0 ) * Hstar );
2273 
2274  double* ysA = ys.data();
2275  double* ytA = yt.data();
2276  double* yrA = yr.data();
2277 
2278  // set x:
2279  for ( j = Nf - 1; j >= 0; j-- )
2280  x .append( ysA[ j ] );
2281 
2282  for ( j = 1; j < Nm; j++ )
2283  x .append( ytA[ j ] );
2284 
2285  for ( j = Js; j >= 0; j-- )
2286  x .append( yrA[ j ] );
2287 
2288  // Smooth out
2289  xA = x.data();
2290  int jj = Nf + Nm;
2291  xA[ jj ] = ( xA[ jj - 1 ] + xA[ jj + 1 ] ) / 2.0;
2292  xA[ jj + 1 ] = ( xA[ jj ] + xA[ jj + 2 ] ) / 2.0;
2293  } // if
2294 }
2295 
2297 //
2298 // mesh_gen_RefL: refine mesh near meniscus (for s>0) or near bottom (for s<0)
2299 // to be used for the acceleration stage
2300 //
2301 // parameters: N0 = number of elements near meniscus (or bottom) to be refined
2302 // M0 = number of elements to be used for the refined region
2303 //
2305 
2306 void US_Astfem_RSA::mesh_gen_RefL( int N0, int M0 )
2307 {
2308  const double PIhalf = M_PI / 2.0;
2309  QVector< double > zz; // temporary array for adaptive grids
2310  double* zA;
2311 
2312  zz.clear();
2313  zz.reserve( x.size() );
2314  xA = x.data();
2315 
2316  if ( US_AstfemMath::minval( af_params.s ) > 0 ) // All species with s>0
2317  {
2318  // Refine around the meniscus for acceleration
2319  for ( int j = 0; j < M0; j++ )
2320  {
2321  double tmp = (double) j / (double)M0;
2322  tmp = 1.0 - cos( tmp * PIhalf );
2323  zz .append( xA[ 0 ] * ( 1.0 - tmp ) + xA[ N0 ] * tmp );
2324  //double tmp = 1.0 - cos( ( (double)j / (double)M0 ) * PIhalf );
2325  //zz << ( xA[ 0 ] * tmp + xA[ N0 ] * ( 1.0 - tmp ) );
2326  }
2327 
2328  for ( int j = N0; j < x.size(); j++ )
2329  //for ( int j = M0; j < x.size(); j++ )
2330  zz .append( xA[ j ] );
2331 
2332 #if 0
2333 double xAval = xA[ x.size() - 1 ];
2334 double xAinc = xAval - xA[ x.size() - 2 ];
2335 zz << xAval + xAinc;
2336 zz << xAval + xAinc * 2.0;
2337 #endif
2338  x.clear();
2339  x.reserve( zz.size() );
2340  zA = zz.data();
2341 
2342  for ( int j = 0; j < zz.size(); j++ )
2343  x .append( zA[ j ] );
2344  }
2345  else if ( US_AstfemMath::maxval( af_params.s ) < 0 ) // All species with s<0
2346  {
2347  for ( int j = 0; j < x.size() - N0; j++ )
2348  zz .append( xA[ j ] );
2349 
2350  // Refine around the bottom for acceleration
2351  int kk = x.size() - 1;
2352  double x1 = xA[ kk - N0 ];
2353  double x2 = xA[ kk ];
2354  double tinc = PIhalf / (double)M0;
2355  double tval = 0.0;
2356 
2357  for ( int j = 1; j <= M0; j++ )
2358  {
2359  tval += tinc;
2360  double tmp = sin( tval );
2361  zz .append( x1 * ( 1.0 - tmp ) + x2 * tmp );
2362  }
2363 
2364  x.clear();
2365  x.reserve( zz.size() );
2366  zA = zz.data();
2367 
2368  for ( int j = 0; j < zz.size(); j++ )
2369  x .append( zA[ j ] );
2370  }
2371  else // Sedimentation and floating mixed up
2372  DbgErr() << "No refinement at ends since sedimentation "
2373  "and floating mixed ...\n" ;
2374 
2375  Nx = x.size();
2376  xA = x.data();
2377 int mm=x.size()/2;
2378 int ee=x.size()-1;
2379 DbgLv(2) << "RSA:mgR: Nx xA sme" << Nx << x[0] << x[1] << x[2]
2380  << x[mm-1] << x[mm] << x[mm+1] << x[mm+2] << x[ee-2] << x[ee-1] << x[ee];
2381 }
2382 
2383 // Compute the coefficient matrices based on fixed mesh
2384 
2385 void US_Astfem_RSA::ComputeCoefMatrixFixedMesh(
2386  double D, double sw2, double** CA, double** CB )
2387 {
2388  if ( Nx != x.size() || Nx < 1 )
2389  DbgErr() << "***FixedMesh ERROR*** Nx x.size" << Nx << x.size()
2390  << " params.s[0] D sw2" << af_params.s[0] << D << sw2;
2391 
2392 #ifdef NO_DB
2393  static int Nsave = 0;
2394  static double*** Stif = NULL;
2395  xA = x.data();
2396 
2397  if ( Nx > Nsave )
2398  {
2399  if ( Nsave > 0 )
2400  US_AstfemMath::clear_3d( Nsave, 4, Stif );
2401 
2402  Nsave = Nx + 200;
2403  US_AstfemMath::initialize_3d( Nsave, 4, 4, &Stif );
2404 //qDebug() << "FixedMesh Nsave" << Nsave;
2405  }
2406 #else
2407  double*** Stif = NULL;
2408  US_AstfemMath::initialize_3d( Nx, 4, 4, &Stif );
2409 #endif
2410 
2411  double xd[ 4 ][ 2 ]; // coord for vertices of quad elem
2412 
2413  for ( int k = 0; k < Nx - 1; k++ )
2414  { // loop for all elem
2415  xd[ 0 ][ 0 ] = xA[ k ];
2416  xd[ 0 ][ 1 ] = 0.0;
2417  xd[ 1 ][ 0 ] = xA[ k + 1 ];
2418  xd[ 1 ][ 1 ] = 0.0;
2419  xd[ 2 ][ 0 ] = xA[ k + 1 ];
2420  xd[ 2 ][ 1 ] = af_params.dt;
2421  xd[ 3 ][ 0 ] = xA[ k ];
2422  xd[ 3 ][ 1 ] = af_params.dt;
2423 
2424  stfb0.CompLocalStif( 4, xd, D, sw2, Stif[ k ] );
2425  }
2426 
2427  // Assemble coefficient matrices
2428  // elem[ 0 ]; i=0
2429  int k = 0;
2430  int m = 0;
2431  CA[ 1 ][ k ] = Stif[ k ][ 3 ][ 0 ] + Stif[ k ][ 3 ][ 3 ]; // j=3;
2432  CA[ 2 ][ k ] = Stif[ k ][ 2 ][ 0 ] + Stif[ k ][ 2 ][ 3 ]; // j=2;
2433  CB[ 1 ][ k ] = Stif[ k ][ 0 ][ 0 ] + Stif[ k ][ 0 ][ 3 ]; // j=0;
2434  CB[ 2 ][ k ] = Stif[ k ][ 1 ][ 0 ] + Stif[ k ][ 1 ][ 3 ]; // j=1;
2435 
2436  for ( k = 1, m = 0; k < Nx - 1; k++, m++ )
2437  { // loop for all elem
2438  // elem k-1: i=1,2
2439  CA[ 0 ][ k ] = Stif[ m ][ 3 ][ 1 ] + Stif[ m ][ 3 ][ 2 ]; // j=3;
2440  CA[ 1 ][ k ] = Stif[ m ][ 2 ][ 1 ] + Stif[ m ][ 2 ][ 2 ]; // j=2;
2441  CB[ 0 ][ k ] = Stif[ m ][ 0 ][ 1 ] + Stif[ m ][ 0 ][ 2 ]; // j=0;
2442  CB[ 1 ][ k ] = Stif[ m ][ 1 ][ 1 ] + Stif[ m ][ 1 ][ 2 ]; // j=1;
2443 
2444  // elem k: i=0,3
2445  CA[ 1 ][ k ] += Stif[ k ][ 3 ][ 0 ] + Stif[ k ][ 3 ][ 3 ]; // j=3;
2446  CA[ 2 ][ k ] = Stif[ k ][ 2 ][ 0 ] + Stif[ k ][ 2 ][ 3 ]; // j=2;
2447  CB[ 1 ][ k ] += Stif[ k ][ 0 ][ 0 ] + Stif[ k ][ 0 ][ 3 ]; // j=0;
2448  CB[ 2 ][ k ] = Stif[ k ][ 1 ][ 0 ] + Stif[ k ][ 1 ][ 3 ]; // j=1;
2449  }
2450 
2451  // elem[ Nx-2 ]; i=1,2
2452  k = Nx - 1;
2453  m = k - 1;
2454  CA[ 0 ][ k ] = Stif[ m ][ 3 ][ 1 ] + Stif[ m ][ 3 ][ 2 ]; // j=3;
2455  CA[ 1 ][ k ] = Stif[ m ][ 2 ][ 1 ] + Stif[ m ][ 2 ][ 2 ]; // j=2;
2456  CB[ 0 ][ k ] = Stif[ m ][ 0 ][ 1 ] + Stif[ m ][ 0 ][ 2 ]; // j=0;
2457  CB[ 1 ][ k ] = Stif[ m ][ 1 ][ 1 ] + Stif[ m ][ 1 ][ 2 ]; // j=1;
2458 
2459 #ifndef NO_DB
2460  US_AstfemMath::clear_3d( Nx, 4, Stif );
2461 #endif
2462 int mm=Nx/2;
2463 DbgLv(2) << "RSA:CCMFM: CA0 sme" << CA[0][0] << CA[0][1] << CA[0][2]
2464  << CA[0][mm-1] << CA[0][mm] << CA[0][mm+1] << CA[0][mm+2]
2465  << CA[0][Nx-3] << CA[0][Nx-2] << CA[0][Nx-1];
2466 DbgLv(2) << "RSA:CCMFM: CA1 sme" << CA[1][0] << CA[1][1] << CA[1][2]
2467  << CA[1][mm-1] << CA[1][mm] << CA[1][mm+1] << CA[1][mm+2]
2468  << CA[1][Nx-3] << CA[1][Nx-2] << CA[1][Nx-1];
2469 DbgLv(2) << "RSA:CCMFM: CB0 sme" << CB[0][0] << CB[0][1] << CB[0][2]
2470  << CB[0][mm-1] << CB[0][mm] << CB[0][mm+1] << CB[0][mm+2]
2471  << CB[0][Nx-3] << CB[0][Nx-2] << CB[0][Nx-1];
2472 DbgLv(2) << "RSA:CCMFM: CB1 sme" << CB[1][0] << CB[1][1] << CB[1][2]
2473  << CB[1][mm-1] << CB[1][mm] << CB[1][mm+1] << CB[1][mm+2]
2474  << CB[1][Nx-3] << CB[1][Nx-2] << CB[1][Nx-1];
2475 }
2476 
2477 void US_Astfem_RSA::ComputeCoefMatrixMovingMeshR(
2478  double D, double sw2, double** CA, double** CB )
2479 {
2480  // Compute local stiffness matrices
2481  double xd[ 4 ][ 2 ]; // coord for verices of quad elem
2482  xA = x.data();
2483 
2484 #ifdef NO_DB
2485  static int Nsave = 0;
2486  static double*** Stif = NULL;
2487 
2488  if ( Nx > Nsave )
2489  {
2490  if ( Nsave > 0 )
2491  US_AstfemMath::clear_3d( Nsave, 4, Stif );
2492 
2493  Nsave = Nx + 200;
2494  US_AstfemMath::initialize_3d( Nsave, 4, 4, &Stif );
2495 //qDebug() << "MovMeshR Nsave" << Nsave;
2496  }
2497 #else
2498  double*** Stif = NULL;
2499  US_AstfemMath::initialize_3d( Nx, 4, 4, &Stif );
2500 #endif
2501 
2502  // elem[0]: triangle
2503  xd[ 0 ][ 0 ] = xA[ 0 ]; xd[ 0 ][ 1 ] = 0.;
2504  xd[ 1 ][ 0 ] = xA[ 1 ]; xd[ 1 ][ 1 ] = af_params.dt;
2505  xd[ 2 ][ 0 ] = xA[ 0 ]; xd[ 2 ][ 1 ] = af_params.dt;
2506  stfb0.CompLocalStif( 3, xd, D, sw2, Stif[ 0 ] );
2507 
2508  // elem[k]: k=1..(Nx-2), quadrilateral
2509  for ( int k = 1; k < Nx - 1; k++ ) // loop for all elem
2510  {
2511  xd[ 0 ][ 0 ] = xA[ k - 1 ]; xd[ 0 ][ 1 ] = 0.0;
2512  xd[ 1 ][ 0 ] = xA[ k ]; xd[ 1 ][ 1 ] = 0.0;
2513  xd[ 2 ][ 0 ] = xA[ k + 1 ]; xd[ 2 ][ 1 ] = af_params.dt;
2514  xd[ 3 ][ 0 ] = xA[ k ]; xd[ 3 ][ 1 ] = af_params.dt;
2515  stfb0.CompLocalStif( 4, xd, D, sw2, Stif[ k ] );
2516  }
2517 
2518  // elem[Nx-1]: triangle
2519  xd[ 0 ][ 0 ] = xA[ Nx - 2 ]; xd[ 0 ][ 1 ] = 0.0;
2520  xd[ 1 ][ 0 ] = xA[ Nx - 1 ]; xd[ 1 ][ 1 ] = 0.0;
2521  xd[ 2 ][ 0 ] = xA[ Nx - 1 ]; xd[ 2 ][ 1 ] = af_params.dt;
2522  stfb0.CompLocalStif( 3, xd, D, sw2, Stif[ Nx - 1 ] );
2523 
2524  // assembly coefficient matrices
2525 
2526  int k = 0;
2527  int mm = 0;
2528  CA[ 1 ][ k ] = Stif[ k ][ 2 ][ 2 ] ;
2529  CA[ 2 ][ k ] = Stif[ k ][ 1 ][ 2 ] ;
2530  CB[ 2 ][ k ] = Stif[ k ][ 0 ][ 2 ] ;
2531 
2532  k = 1;
2533  CA[ 0 ][ k ] = Stif[ mm ][ 2 ][ 0 ] + Stif[ mm ][ 2 ][ 1 ];
2534  CA[ 1 ][ k ] = Stif[ mm ][ 1 ][ 0 ] + Stif[ mm ][ 1 ][ 1 ];
2535  CA[ 1 ][ k ] += Stif[ k ][ 3 ][ 0 ] + Stif[ k ][ 3 ][ 3 ];
2536  CA[ 2 ][ k ] = Stif[ k ][ 2 ][ 0 ] + Stif[ k ][ 2 ][ 3 ] ;
2537 
2538  CB[ 1 ][ k ] = Stif[ mm ][ 0 ][ 0 ] + Stif[ mm ][ 0 ][ 1 ]; // j=0;
2539  CB[ 1 ][ k ] += Stif[ k ][ 0 ][ 0 ] + Stif[ k ][ 0 ][ 3 ]; // j=0;
2540  CB[ 2 ][ k ] = Stif[ k ][ 1 ][ 0 ] + Stif[ k ][ 1 ][ 3 ]; // j=1;
2541 
2542  for( k = 2; k < Nx - 1; k++ )
2543  { // loop for all elem
2544  // elem k-1: i=1,2
2545  mm = k - 1;
2546  CA[ 0 ][ k ] = Stif[ mm ][ 3 ][ 1 ] + Stif[ mm ][ 3 ][ 2 ]; // j=3;
2547  CA[ 1 ][ k ] = Stif[ mm ][ 2 ][ 1 ] + Stif[ mm ][ 2 ][ 2 ]; // j=2;
2548  CB[ 0 ][ k ] = Stif[ mm ][ 0 ][ 1 ] + Stif[ mm ][ 0 ][ 2 ]; // j=0;
2549  CB[ 1 ][ k ] = Stif[ mm ][ 1 ][ 1 ] + Stif[ mm ][ 1 ][ 2 ]; // j=1;
2550 
2551  // elem k: i=0,3
2552  CA[ 1 ][ k ] += Stif[ k ][ 3 ][ 0 ] + Stif[ k ][ 3 ][ 3 ]; // j=3;
2553  CA[ 2 ][ k ] = Stif[ k ][ 2 ][ 0 ] + Stif[ k ][ 2 ][ 3 ]; // j=2;
2554  CB[ 1 ][ k ] += Stif[ k ][ 0 ][ 0 ] + Stif[ k ][ 0 ][ 3 ]; // j=0;
2555  CB[ 2 ][ k ] = Stif[ k ][ 1 ][ 0 ] + Stif[ k ][ 1 ][ 3 ]; // j=1;
2556  }
2557 
2558  k = Nx - 1;
2559  mm = k - 1;
2560  // elem[k-1]: quadrilateral
2561  CA[ 0 ][ k ] = Stif[ mm ][3][1] + Stif[ mm ][3][2]; // j=3;
2562  CA[ 1 ][ k ] = Stif[ mm ][2][1] + Stif[ mm ][2][2]; // j=2;
2563  CB[ 0 ][ k ] = Stif[ mm ][0][1] + Stif[ mm ][0][2]; // j=0;
2564  CB[ 1 ][ k ] = Stif[ mm ][1][1] + Stif[ mm ][1][2]; // j=1;
2565 
2566  // elem[k]: triangle
2567  CA[ 1 ][ k ] += Stif[ k ][2][0] + Stif[ k ][2][1] + Stif[ k ][2][2];
2568  CB[ 1 ][ k ] += Stif[ k ][0][0] + Stif[ k ][0][1] + Stif[ k ][0][2];
2569  CB[ 2 ][ k ] = Stif[ k ][1][0] + Stif[ k ][1][1] + Stif[ k ][1][2];
2570 
2571 #ifndef NO_DB
2572  US_AstfemMath::clear_3d( Nx, 4, Stif );
2573 #endif
2574 }
2575 
2576 void US_Astfem_RSA::ComputeCoefMatrixMovingMeshL(
2577  double D, double sw2, double** CA, double** CB )
2578 {
2579  // compute local stiffness matrices
2580  double xd[4][2]; // coord for verices of quad elem
2581  xA = x.data();
2582 
2583 #ifdef NO_DB
2584  static int Nsave = 0;
2585  static double*** Stif = NULL;
2586 
2587  if ( Nx > Nsave )
2588  {
2589  if ( Nsave > 0 )
2590  US_AstfemMath::clear_3d( Nsave, 4, Stif );
2591 
2592  Nsave = Nx + 200;
2593  US_AstfemMath::initialize_3d( Nsave, 4, 4, &Stif );
2594  }
2595 #else
2596  double*** Stif = NULL;
2597  US_AstfemMath::initialize_3d( Nx, 4, 4, &Stif );
2598 #endif
2599 
2600  // elem[0]: triangle
2601  xd[0][0] = xA[0];
2602  xd[0][1] = 0.0;
2603  xd[1][0] = xA[1]; xd[1][1] = 0.0;
2604  xd[2][0] = xA[0]; xd[2][1] = af_params.dt;
2605  stfb0.CompLocalStif( 3, xd, D, sw2, Stif[ 0 ] );
2606 
2607  // elem[k]: k=1..(Nx-2), quadrilateral
2608  for ( int k = 1; k < Nx - 1; k++ )
2609  { // loop for all elem
2610  xd[0][0] = xA[k ]; xd[0][1] = 0.0;
2611  xd[1][0] = xA[k+1]; xd[1][1] = 0.0;
2612  xd[2][0] = xA[k ]; xd[2][1] = af_params.dt;
2613  xd[3][0] = xA[k-1]; xd[3][1] = af_params.dt;
2614  stfb0.CompLocalStif( 4, xd, D, sw2, Stif[ k ] );
2615  }
2616 
2617  // elem[Nx-1]: triangle
2618  xd[0][0] = xA[Nx-1]; xd[0][1] = 0.0;
2619  xd[1][0] = xA[Nx-1]; xd[1][1] = af_params.dt;
2620  xd[2][0] = xA[Nx-2]; xd[2][1] = af_params.dt;
2621  stfb0.CompLocalStif( 3, xd, D, sw2, Stif[ Nx - 1 ] );
2622 
2623  // assembly coefficient matrices
2624 
2625  int k = 0;
2626  CA[1][0] = Stif[0][2][0] + Stif[0][2][1] + Stif[0][2][2];
2627  CB[0][0] = Stif[0][0][0] + Stif[0][0][1] + Stif[0][0][2] ;
2628  CB[1][0] = Stif[0][1][0] + Stif[0][1][1] + Stif[0][1][2] ;
2629 
2630  CA[1][0]+= Stif[1][3][0] + Stif[1][3][3] ;
2631  CA[2][0] = Stif[1][2][0] + Stif[1][2][3] ;
2632  CB[1][0]+= Stif[1][0][0] + Stif[1][0][3] ;
2633  CB[2][0] = Stif[1][1][0] + Stif[1][1][3] ;
2634 
2635  for ( int k = 1; k < Nx - 2; k++ )
2636  { // loop for all elem
2637  // elem k:
2638  CA[0][k] = Stif[k ][3][1] + Stif[k ][3][2]; // j=3;
2639  CA[1][k] = Stif[k ][2][1] + Stif[k ][2][2]; // j=2;
2640  CB[0][k] = Stif[k ][0][1] + Stif[k ][0][2]; // j=0;
2641  CB[1][k] = Stif[k ][1][1] + Stif[k ][1][2]; // j=1;
2642 
2643  // elem k+1:
2644  CA[1][k] += Stif[k+1][3][0] + Stif[k+1][3][3]; // j=3;
2645  CA[2][k] = Stif[k+1][2][0] + Stif[k+1][2][3]; // j=2;
2646  CB[1][k] += Stif[k+1][0][0] + Stif[k+1][0][3]; // j=0;
2647  CB[2][k] = Stif[k+1][1][0] + Stif[k+1][1][3]; // j=1;
2648  }
2649 
2650  k = Nx - 2;
2651  // elem k:
2652  CA[0][k] = Stif[k ][3][1] + Stif[k ][3][2]; // j=3;
2653  CA[1][k] = Stif[k ][2][1] + Stif[k ][2][2]; // j=2;
2654  CB[0][k] = Stif[k ][0][1] + Stif[k ][0][2]; // j=0;
2655  CB[1][k] = Stif[k ][1][1] + Stif[k ][1][2]; // j=1;
2656 
2657  // elem k+1: (triangle)
2658  CA[1][k] += Stif[k+1][2][0] + Stif[k+1][2][2]; // j=3;
2659  CA[2][k] = Stif[k+1][1][0] + Stif[k+1][1][2]; // j=2;
2660  CB[1][k] += Stif[k+1][0][0] + Stif[k+1][0][2]; // j=0;
2661 
2662  k = Nx - 1;
2663  // elem[k]: triangle
2664  CA[0][k] = Stif[k ][2][1] ;
2665  CA[1][k] = Stif[k ][1][1] ;
2666  CB[0][k] = Stif[k ][0][1] ;
2667 
2668 #ifndef NO_DB
2669  US_AstfemMath::clear_3d( Nx, 4, Stif );
2670 #endif
2671 }
2672 
2673 // Given total concentration of a group of components involved,
2674 // find the concentration of each component by equilibrium condition
2675 void US_Astfem_RSA::decompose( US_AstfemMath::MfemInitial* C0 )
2676 {
2677  int num_comp = af_params.role.size();
2678 
2679  // Note: all components must be defined on the same radial grids
2680  int Npts = C0[ 0 ].radius.size();
2681 //DbgLv(2) << "RSA: decompose() num_comp Npts" << num_comp << Npts;
2682 
2683  // Special case: self-association n A <--> An
2684  if ( num_comp == 2 ) // Only 2 components and one association rule
2685  {
2686  int st0 = af_params.association[ 0 ].stoichs[ 0 ];
2687  int st1 = af_params.association[ 0 ].stoichs[ 1 ];
2688  double k_d = af_params.association[ 0 ].k_d;
2689  double k_assoc = ( k_d != 0.0 ) ? ( 1.0 / k_d ) : 0.0;
2690  // Extinction coefficient for the monomer, corrected for pathlength
2691  double ext_M = af_params.kext[ 0 ];
2692  // K_association converted from signal to molar units
2693  k_assoc /= ext_M;
2694 DbgLv(2) << "RSA: decompose() ext_M" << ext_M << "k_assoc" << k_assoc;
2695 #ifndef NO_DB
2696  emit current_component( -Npts );
2697 #endif
2698 
2699  for ( int j = 0; j < Npts; j++ )
2700  {
2701  double c1 = 0.0;
2702  double ct = C0[ 0 ].concentration[ j ] + C0[ 1 ].concentration[ j ] ;
2703 //DbgLv(2) << "RSA: j st0 st1" << j << st0 << st1;
2704 
2705  if ( st0 == 2 && st1 == -1 ) // mono <--> dimer
2706  {
2707  //double ext_Msq = ext_M * ext_M;
2708  //c1 = ( sqrt( ext_Msq + 4.0 * k_assoc * ct * ext_Msq) - ext_M ) / ( 2.0 * k_assoc * ext_Msq);
2709  c1 = ( sqrt( 1.0 + 4.0 * k_assoc * ct ) - 1.0 ) / ( 2.0 * k_assoc );
2710  }
2711 
2712  else if ( st0 == 3 && st1 == -1 ) // mono <--> trimer
2713  c1 = US_AstfemMath::cube_root( -ct / k_assoc, 1.0 / k_assoc, 0.0 );
2714 
2715  else if ( st0 > 3 && st1 == -1 ) // mono <--> n-mer
2716  c1 = US_AstfemMath::find_C1_mono_Nmer( st0, k_assoc, ct );
2717 
2718  else
2719  {
2720  DbgErr() << "Warning: invalid stoichiometry in decompose()"
2721  << " st0 st1 c1" << st0 << st1 << c1;
2722  return;
2723  }
2724 
2725  double c2 = k_assoc * pow( c1, (double)st0 );
2726 if(j<5||(j+6)>Npts)
2727 DbgLv(2) << "RSA: decompose() j" << j << "c1 c2 ct Ka" << c1 << c2 << ct
2728  << k_assoc << "p_role0_st0" << af_params.role[0].stoichs[0];
2729 
2730  if ( af_params.role[ 0 ].stoichs[ 0 ] > 0 ) // c1=reactant
2731  {
2732  C0[ 0 ].concentration[ j ] = c1 ;
2733  C0[ 1 ].concentration[ j ] = c2 ;
2734  }
2735  else // c1=product
2736  {
2737  C0[ 0 ].concentration[ j ] = c2 ;
2738  C0[ 1 ].concentration[ j ] = c1 ;
2739  }
2740 #ifndef NO_DB
2741  emit current_component( j + 1 );
2742 #endif
2743  qApp->processEvents();
2744  if ( stopFlag ) break;
2745  }
2746 DbgLv(2) << "RSA: decompose NCOMP=2 return";
2747  return;
2748  }
2749 
2750  // General cases
2751  double** C1;
2752  double** C2; // Arrays for all components at all radius position
2753 
2754  US_AstfemMath::initialize_2d( num_comp, Npts, &C1 );
2755  US_AstfemMath::initialize_2d( num_comp, Npts, &C2 );
2756 
2757  for( int i = 0; i < num_comp; i++ )
2758  {
2759  for( int j = 0; j < Npts; j++ )
2760  {
2761  C1[ i ][ j ] = C0[ i ].concentration[ j ];
2762  C2[ i ][ j ] = C1[ i ][ j ];
2763  }
2764  }
2765 
2766  // Estimate max time to reach equilibrium and suitable step size:
2767  // using e^(-k_min * Nx * dt ) < 1.e-7
2768 
2769  double k_min = 1.0e12;
2770 
2771  // Get minimum k
2772  for ( int i = 0; i < af_params.association.size(); i++ )
2773  {
2774  k_min = qMin( k_min, af_params.association[ i ].k_off );
2775  }
2776 
2777  k_min = qMin( k_min, 1.0e-12 );
2778 
2779  // Max number of time steps to get to equilibrium
2780  const int time_max = 100;
2781  double timeStepSize = - log( 1.0e-7 ) / ( k_min * time_max );
2782 DbgLv(2) << "RSA: decompose k_min time_max timeStepSize"
2783  << k_min << time_max << timeStepSize;
2784 
2785  // time loop
2786 #ifndef NO_DB
2787  emit calc_start( time_max );
2788  emit current_component( -time_max );
2789 #endif
2790 
2791  for ( int ti = 0; ti < time_max; ti++ )
2792  {
2793 #ifndef NO_DB
2794  if ( show_movie && (ti%8) == 0 )
2795  {
2796 //DbgLv(2) << "AR: calc_progr ti" << ti;
2797  emit calc_progress( ti );
2798  qApp->processEvents();
2799  //US_Sleep::msleep( 10 );
2800  }
2801 #endif
2802 
2803  ReactionOneStep_Euler_imp( Npts, C1, timeStepSize );
2804 
2805  if ( stopFlag ) break;
2806 
2807  double diff = 0.0;
2808  double ct = 0.0;
2809 
2810  for ( int i = 0; i < num_comp; i++ )
2811  {
2812  for ( int j = 0; j < Npts; j++ )
2813  {
2814  diff += qAbs( C2[ i ][ j ] - C1[ i ][ j ] );
2815  ct += qAbs( C1[ i ][ j ] );
2816  C2[ i ][ j ] = C1[ i ][ j ];
2817  }
2818  }
2819 
2820 #ifndef NO_DB
2821  emit current_component( ti + 1 );
2822  qApp->processEvents();
2823 #endif
2824 
2825  if ( diff < 1.0e-5 * ct )
2826  {
2827 #ifndef NO_DB
2828  int step = ti + 1;
2829  emit current_component( -step );
2830  emit current_component( step );
2831 #endif
2832  break;
2833  }
2834  } // end time steps
2835 
2836 #ifndef NO_DB
2837  if ( show_movie )
2838  {
2839  emit calc_done();
2840  qApp->processEvents();
2841  }
2842 #endif
2843 
2844  for ( int i = 0; i < num_comp; i++ )
2845  {
2846  for ( int j = 0; j < Npts; j++ )
2847  C0[ i ].concentration[ j ] = C1[ i ][ j ];
2848  }
2849 
2850  US_AstfemMath::clear_2d( num_comp, C1 );
2851  US_AstfemMath::clear_2d( num_comp, C2 );
2852 }
2853 
2855 //
2856 // ReactionOneStep_Euler_imp: implicit Mid-point Euler
2857 //
2859 void US_Astfem_RSA::ReactionOneStep_Euler_imp(
2860  int Npts, double** C1, double timeStep )
2861 {
2862  int num_comp = af_params.role.size();
2863 DbgLv(2) << "RSA:Eul: Npts timeStep" << Npts << timeStep
2864  << "num_comp" << num_comp;
2865 
2866  // Special case: self-association n A <--> An
2867  if ( num_comp == 2 ) // only 2 components and one association rule
2868  {
2869  double uhat;
2870 
2871  // Current rule used
2872  int rule = rg[ af_params.rg_index ].association[ 0 ];
2873  int st0 = system.associations[ rule ].stoichs[ 0 ];
2874  int st1 = system.associations[ rule ].stoichs[ 1 ];
2875  double extn = af_params.kext[ system.associations[rule].rcomps[0] ];
2876  double k_d = system.associations[ rule ].k_d;
2877  double k_assoc = ( k_d == 0.0 ) ? 0.0 : ( 1.0 / ( k_d * extn ) );
2878  double k_off = system.associations[ rule ].k_off;
2879 DbgLv(2) << "RSA:Eul: rule" << rule << "st0 st1 k_assoc k_off"
2880  << st0 << st1 << k_assoc << k_off;
2881 
2882  for ( int j = 0; j < Npts; j++ )
2883  {
2884  double ct = C1[ 0 ][ j ] + C1[ 1 ][ j ];
2885  double dva = timeStep * k_off * k_assoc;
2886  double dvb = timeStep * k_off + 2.;
2887  double dvc = timeStep * k_off * ct + 2.0 * C1[ 0 ][ j ];
2888 if((j+9)>Npts)
2889 DbgLv(2) << "RSA:Eul: j ct" << j << ct << "dva dvb dvc" << dva << dvb << dvc;
2890 
2891  if ( st0 == 2 && st1 == (-1) ) // mono <--> dimer
2892  {
2893  uhat = 2. * dvc / ( dvb + sqrt( dvb * dvb + 4. * dva * dvc ) );
2894 if((j+9)>Npts)
2895 DbgLv(2) << "RSA:Eul: mono-dimer uhat" << uhat
2896  << "C1[0][j] C1[1][j]" << C1[0][j] << C1[1][j];
2897  }
2898 
2899  else if ( st0 == 3 && st1 == (-1) ) // mono <--> trimer
2900  {
2901  uhat = US_AstfemMath::cube_root( -dvc / dva, dvb / dva, 0.0 );
2902  }
2903 
2904  else if ( st0 > 3 && st1 == (-1) ) // mono <--> n-mer
2905  {
2906  uhat = US_AstfemMath::find_C1_mono_Nmer( st0, dva / dvb, dvc / dvb);
2907  }
2908 
2909  else
2910  {
2911  DbgErr() << "Warning: invalid stoichiometry in decompose()";
2912  return;
2913  }
2914 
2915  if ( af_params.role[ 0 ].stoichs[ 0 ] > 0 ) // c1=reactant
2916  {
2917  C1[ 0 ][ j ] = 2. * uhat - C1[ 0 ][ j ];
2918  C1[ 1 ][ j ] = ct - C1[ 0 ][ j ];
2919 //DbgLv(2) << "RSA:Eul: c1=react c1[0][j] c1[1][j]" << C1[0][j] << C1[1][j];
2920  }
2921 
2922  else
2923  { // c1=product
2924  C1[ 1 ][ j ] = 2. * uhat - C1[ 1 ][ j ];
2925  C1[ 0 ][ j ] = ct - C1[ 1 ][ j ];
2926 //DbgLv(2) << "RSA:Eul: c1=prod c1[0][j] c1[1][j]" << C1[0][j] << C1[1][j];
2927  }
2928  }
2929  return;
2930  }
2931 
2932  // General cases
2933  const int iter_max = 20; // maximum number of Newton iteration allowed
2934 
2935  double** A;
2936 
2937  QVector< double > y0Vec( num_comp );
2938  QVector< double > dnVec( num_comp );
2939  QVector< double > bbVec( num_comp );
2940  QVector< double > y0rVc( num_comp );
2941  QVector< double > y1rVc( num_comp );
2942  y0rVc.fill( 0.0, num_comp );
2943  y1rVc.fill( 0.0, num_comp );
2944  double* y0 = y0Vec.data();
2945  double* delta_n = dnVec.data();
2946  double* b = bbVec.data();
2947  double* y0_ref = y0rVc.data();
2948  double* y1_ref = y1rVc.data();
2949 
2950  US_AstfemMath::initialize_2d( num_comp, num_comp, &A );
2951 
2952  for ( int j = 0; j < Npts; j++ )
2953  {
2954  double ct = 0.0;
2955  double diff_ref = 0.0;
2956 
2957  for ( int i = 0; i < num_comp; i++ )
2958  {
2959  y0[ i ] = C1[ i ][ j ];
2960  delta_n[ i ] = 0.0;
2961  ct += qAbs( y0[ i ] );
2962  diff_ref += qAbs( y0_ref[ i ] - y0[ i ] );
2963  }
2964 
2965  if ( diff_ref < ( ct * 1.e-7 ) || diff_ref < 1.e-9 )
2966  {
2967  for ( int i = 0; i < num_comp; i++ )
2968  {
2969  C1[ i ][ j ] = y1_ref[ i ];
2970  }
2971 
2972  continue;
2973  }
2974 
2975  for ( int iter = 0; iter < iter_max; iter++ ) // Newton iteration
2976  {
2977  double diff = 0.0;
2978 
2979  for ( int i = 0; i < num_comp; i++ )
2980  {
2981  y0[ i ] = C1[ i ][ j ] + delta_n[ i ];
2982  }
2983 
2984 //qDebug() << "RSA:Eul: j" << j << "iter" << iter << " y0[0]" << y0[0];
2985  Reaction_dydt( y0, b ); // b=dy/dt
2986 //qDebug() << "RSA:Eul: post-dydt y0[0]" << y0[0] << "b0" << b[0];
2987  Reaction_dfdy( y0, A ); // A=df/dy
2988 //qDebug() << "RSA:Eul: post-dfdy y0[0]" << y0[0] << "A00" << A[0][0]
2989 // << "A[n][n]" << A[num_comp-1][num_comp-1];
2990 
2991  for ( int i = 0; i < num_comp; i++ )
2992  {
2993  for ( int k = 0; k < num_comp; k++ )
2994  A[ i ][ k ] *= ( -timeStep );
2995 
2996  A[ i ][ i ] += 2.0;
2997  //b[ i ] = - ( 2 * delta_n[ i ] - timeStep * b[ i ] );
2998  b[ i ] = b[ i ] * timeStep - delta_n[ i ] * 2.0;
2999  }
3000 
3001  if ( US_AstfemMath::GaussElim( num_comp, A, b ) == -1 )
3002  {
3003  DbgErr() << "Matrix singular in Reaction_Euler_imp: model 12";
3004  break;
3005  }
3006  else
3007  {
3008  diff = 0.0;
3009 
3010  for ( int i = 0; i < num_comp; i++ )
3011  {
3012  delta_n[ i ] += b[ i ];
3013  diff += qAbs( delta_n[ i ] );
3014  }
3015  }
3016 
3017  if ( diff < ( 1.0e-7 * ct ) ) break;
3018  } // End of Newton iteration;
3019 
3020  for ( int i = 0; i < num_comp; i++ )
3021  {
3022  y0_ref[ i ] = C1[ i ][ j ];
3023  C1[ i ][ j ] += delta_n[ i ];
3024  y1_ref[ i ] = C1[ i ][ j ];
3025  }
3026 //qDebug() << "RSA:Eul: j" << j << "ct diff_ref" << ct << diff_ref;
3027 
3028  } // End of j (pts)
3029 //int nn=num_comp-1;
3030 //int mm=Npts-1;
3031 //DbgLv(2) << "RSA: ReaEul: num_comp Npts" << num_comp << Npts;
3032 //DbgLv(2) << "RSA: ReaEul: C1[0][0] C1[N][M]" << C1[0][0] << C1[nn][mm]
3033 // << "delta_n[N]" << delta_n[nn] << "b[N]" << b[nn];
3034 
3035  US_AstfemMath::clear_2d( num_comp, A );
3036 }
3037 
3038 void US_Astfem_RSA::Reaction_dydt( double* y0, double* yt )
3039 {
3040  US_AstfemMath::ReactionGroup* rgp = &rg[ af_params.rg_index ];
3041  int num_comp = rgp->GroupComponent.size();
3042  int num_rule = rgp->association.size();
3043  QVector< double > qqVec( num_rule );
3044  double* Q = qqVec.data();
3045 
3046  for ( int m = 0; m < num_rule; m++ )
3047  {
3048  US_Model::Association* as = &af_params.association[ m ];
3049  double k_off = as->k_off;
3050  double k_assoc = ( as->k_d != 0.0 ) ? ( 1.0 / as->k_d ) : 0.0;
3051  k_assoc /= af_params.kext[ as->rcomps[ 0 ] ];
3052  double k_on = k_assoc * k_off;
3053  double Q_reactant = 1.0;
3054  double Q_product = 1.0;
3055 //DbgLv(2) << "RSA:ReDydt: m k_off k_on" << m << k_off << k_on;
3056 
3057  for ( int n = 0; n < as->rcomps.size(); n++ )
3058  {
3059  // local index of the n-th component in assoc[rule]
3060  int ind_cn = as->rcomps[ n ] ;
3061 
3062  // stoichiometry of n-th component in the rule
3063  int kstoi = as->stoichs[ n ] ;
3064  int react = ( kstoi < 0 ) ? -1 : 1;
3065  double rstoi = (double)qAbs( kstoi );
3066 
3067  // extinction coefficient of n'th component
3068  double extn = af_params.kext[ rgp->GroupComponent[ ind_cn ] ];
3069 //if(m==0)
3070 //DbgLv(2) << "RSA:ReDydt: n ind_cn rstoi react extn" << n << ind_cn << rstoi
3071 // << react << extn;
3072 
3073  if ( react > 0 ) // comp[n] here is reactant
3074  {
3075  Q_reactant *= pow( y0[ ind_cn ] / extn, rstoi );
3076  }
3077  else // comp[n] here is product
3078  {
3079  Q_product *= pow( y0[ ind_cn ] / extn, rstoi );
3080  }
3081  }
3082 
3083  Q[ m ] = k_on * Q_reactant - k_off * Q_product;
3084  }
3085 
3086  for ( int i = 0; i < num_comp; i++ )
3087  {
3088  yt[ i ] = 0.0;
3089  US_AstfemMath::ComponentRole* cr = &af_params.role[ i ];
3090 
3091  for ( int m = 0; m < cr->assocs.size(); m++ )
3092  {
3093  yt[ i ] -= ( (double)cr->stoichs[ m ] * Q[ cr->assocs[ m ] ] );
3094  }
3095 
3096  // convert molar into signal concentration
3097  yt[ i ] *= af_params.kext[ rgp->GroupComponent[ i ] ];
3098  }
3099 DbgLv(2) << "RSA:ReDydt: yt[0] yt[n]" << yt[0] << yt[num_comp-1];
3100 }
3101 
3102 void US_Astfem_RSA::Reaction_dfdy( double* y0, double** dfdy )
3103 {
3104  double** QC;
3105 
3106  US_AstfemMath::ReactionGroup* rgp = &rg[ af_params.rg_index ];
3107  int num_comp = rgp->GroupComponent.size();
3108  int num_rule = rgp->association.size();
3109 
3110  US_AstfemMath::initialize_2d( num_rule, num_comp, &QC );
3111 
3112  for ( int m = 0; m < num_rule; m++ )
3113  {
3114  US_Model::Association* as = &af_params.association[ m ];
3115  double k_off = as->k_off;
3116  double k_assoc = ( as->k_d != 0.0 ) ? ( 1.0 / as->k_d ) : 0.0;
3117  k_assoc /= af_params.kext[ as->rcomps[ 0 ] ];
3118  double k_on = k_assoc * k_off;
3119 
3120  for ( int j = 0; j < num_comp; j++ )
3121  {
3122  double Q_reactant = 1.0;
3123  double Q_product = 1.0;
3124  double deriv_r = 0.0;
3125  double deriv_p = 0.0;
3126 
3127  for( int n = 0; n < as->rcomps.size(); n++ )
3128  {
3129  // Local index of the n-th component in assoc[rule]
3130  int ind_cn = as->rcomps[ n ] ;
3131 
3132  // Stoichiometry of n-th comp in the rule
3133  int kstoi = as->stoichs[ n ] ;
3134  int react = ( kstoi < 0 ) ? -1 : 1;
3135  double rstoi = (double)( kstoi * react );
3136 
3137  // Extinction coefficient of n'th component
3138  double extn = af_params.kext[ rgp->GroupComponent[ ind_cn ] ];
3139  double yext = y0[ ind_cn ] / extn;
3140 
3141  // comp[j] is in the rule
3142  if ( as->rcomps[ n ] == j )
3143  {
3144  if ( react > 0 ) // comp[n] is reactant
3145  deriv_r = rstoi / extn * pow( yext, rstoi - 1.0 );
3146 
3147  else // comp[n] in this rule is product
3148  deriv_p = rstoi / extn * pow( yext, rstoi - 1.0 );
3149  }
3150 
3151  else // comp[j] is not in the rule
3152  {
3153  if ( react > 0 ) // comp[n] is reactant
3154  Q_reactant *= pow( yext, rstoi );
3155 
3156  else // comp[n] in this rule is product
3157  Q_product *= pow( yext, rstoi );
3158  }
3159  }
3160 
3161  QC[ m ][ j ] = k_on * Q_reactant * deriv_r
3162  - k_off * Q_product * deriv_p;
3163  } // C_j
3164  } // m-rule
3165 
3166  for ( int i = 0; i < num_comp; i++ )
3167  {
3168  US_AstfemMath::ComponentRole* cr = &af_params.role[ i ];
3169 
3170  for ( int j = 0; j < num_comp; j++ )
3171  {
3172  dfdy[ i ][ j ] = 0.0;
3173 
3174  for ( int m = 0; m < cr->assocs.size(); m++ )
3175  {
3176  dfdy[ i ][ j ] -= ( (double)cr->stoichs[ m ]
3177  * QC[ cr->assocs[ m ] ][ j ] );
3178  }
3179 
3180  // convert molar into signal concentration
3181  dfdy[ i ][ j ] *= af_params.kext[ rgp->GroupComponent[ i ] ];
3182  }
3183  }
3184 
3185  US_AstfemMath::clear_2d( num_rule, QC );
3186 }
3187 
3188 // *** this is the SNI version of operator scheme
3189 
3190 int US_Astfem_RSA::calculate_ra2( double rpm_start, double rpm_stop,
3191  US_AstfemMath::MfemInitial* C_init, US_AstfemMath::MfemData& simdata,
3192  bool accel )
3193 {
3194  int Mcomp = af_params.s.size();
3195 
3196  US_AstfemMath::MfemScan simscan;
3197 
3198  // Generate the adaptive mesh
3199  QVector< double > nu;
3200  nu.clear();
3201  nu.reserve( Mcomp );
3202 
3203  for ( int i = 0; i < Mcomp; i++ )
3204  {
3205  double sw2 = af_params.s[ i ] * sq( rpm_stop * M_PI / 30 );
3206  nu .append( sw2 / af_params.D[ i ] );
3207  }
3208 
3209  mesh_gen( nu, simparams.meshType );
3210 
3211  simdata.radius.clear();
3212  simdata.scan .clear();
3213  simdata.radius.reserve( Nx );
3214  simdata.scan .reserve( Nx );
3215  xA = x.data();
3216 DbgLv(2) << "RSA:_ra2: Mcomp Nx" << Mcomp << Nx << "x size" << x.size()
3217  << "D0 nu0 D1 nu1" << af_params.D[0] << nu[0] << af_params.D[1] << nu[1];
3218 
3219  bool fixedGrid = ( simparams.gridType == US_SimulationParameters::FIXED );
3220  double meniscus = af_params.current_meniscus;
3221  double bottom = af_params.current_bottom;
3222  int Nt = af_params.time_steps;
3223 
3224  // Refine left hand side (when s_max>0) or
3225  // right hand side (when s<0) for acceleration
3226 
3227  // Used for mesh and dt
3228  double s_max = US_AstfemMath::maxval( af_params.s );
3229  double s_min = US_AstfemMath::minval( af_params.s );
3230 
3231  if ( accel )
3232  {
3233  double xc ;
3234 DbgLv(2) << "RSA:_ra2:(2) Nx" << Nx << "x size" << x.size();
3235 
3236  if ( s_min > 0.0 ) // all sediment towards bottom
3237  {
3238  int j;
3239  double sw2 = s_max * sq( rpm_stop * M_PI / 30. );
3240  xc = meniscus + sw2 * ( Nt * af_params.dt ) / 3.;
3241 
3242  for ( j = 0; j < Nx - 3; j++ )
3243  {
3244  if ( xA[ j ] > xc ) break;
3245  }
3246 
3247  mesh_gen_RefL( j + 1, 4 * j );
3248 DbgLv(2) << "RSA:_ra2: s_min s_max" << s_min << s_max << "xc xAj"
3249  << xc << xA[j] << "j" << j << "N0 M0" << (j+1) << (j*4);
3250  }
3251  else if ( s_max < 0.0 ) // all float towards meniscus
3252  {
3253  // s_min corresponds to fastest component
3254  int j;
3255  double sw2 = s_min * sq( rpm_stop * M_PI / 30. );
3256  xc = bottom + sw2 * ( Nt * af_params.dt ) / 3.;
3257 
3258  for ( j = Nx - 1; j > 1; j-- )
3259  {
3260  if ( xA[ j ] < xc ) break;
3261  }
3262 
3263  mesh_gen_RefL( j + 1, 4 * j );
3264  }
3265  else
3266  {
3267  DbgErr() << "Multicomponent system with sedimentation and "
3268  "floating mixed, use uniform mesh";
3269  }
3270 DbgLv(2) << "RSA:_ra2:(3) Nx" << Nx << "x size" << x.size();
3271  }
3272 DbgLv(2) << "RSA:_ra2:(4) Nx" << Nx << "x size" << x.size();
3273 
3274  for ( int i = 0; i < Nx; i++ )
3275  simdata.radius .append( xA[ i ] );
3276 
3277  // Stiffness matrix on left hand side
3278  // CA[0...Ms-1][4][0...Nx-1]
3279  double*** CA;
3280  double*** CA1;
3281  double*** CA2;
3282 
3283  // Stiffness matrix on right hand side
3284  // CB[0...Ms-1][4][0...Nx-1]
3285  double*** CB;
3286  double*** CB1;
3287  double*** CB2;
3288 
3289  // Initialize the coefficient matrices
3290  US_AstfemMath::initialize_3d( Mcomp, 4, Nx, &CA );
3291  US_AstfemMath::initialize_3d( Mcomp, 4, Nx, &CB );
3292 
3293  if ( accel ) // Acceleration, so use fixed grid
3294  {
3295  US_AstfemMath::initialize_3d( Mcomp, 3, Nx, &CA1 );
3296  US_AstfemMath::initialize_3d( Mcomp, 3, Nx, &CA2 );
3297  US_AstfemMath::initialize_3d( Mcomp, 3, Nx, &CB1 );
3298  US_AstfemMath::initialize_3d( Mcomp, 3, Nx, &CB2 );
3299 
3300  for( int i = 0; i < Mcomp; i++ )
3301  {
3302  double sw2 = 0.0;
3303  ComputeCoefMatrixFixedMesh( af_params.D[ i ], sw2, CA1[ i ], CB1[ i ] );
3304 
3305  sw2 = af_params.s[ i ] * sq( rpm_stop * M_PI / 30 );
3306  ComputeCoefMatrixFixedMesh( af_params.D[ i ], sw2, CA2[ i ], CB2[ i ] );
3307  }
3308 DbgLv(2) << "RSA:_ra2:(5) Nx" << Nx << "x size" << x.size();
3309  }
3310 
3311  else // Constant sedimentation speed
3312  {
3313  if ( fixedGrid )
3314  {
3315  for( int i = 0; i < Mcomp; i++ )
3316  {
3317  double sw2 = af_params.s[ i ] * sq( rpm_stop * M_PI / 30 );
3318  ComputeCoefMatrixFixedMesh( af_params.D[ i ], sw2, CA[ i ], CB[ i ] );
3319  }
3320  }
3321 
3322  else // Moving grid
3323  {
3324  if ( s_min > 0) // All components sedimenting
3325  {
3326  double stop_fact = sq( rpm_stop * M_PI / 30.0 );
3327  double sqb = sq( af_params.current_bottom );
3328 
3329  for ( int i = 0; i < Mcomp; i++ )
3330  {
3331  double sw2 = af_params.s[ i ] * stop_fact;
3332  double sw2D = 0.5 * sw2 / af_params.D[ i ];
3333  double s_rat = af_params.s[ i ] / s_max;
3334 DbgLv(2) << "RSA: smin>0:GlStf: i" << i << "Si Di sw2 sw2D sqb s_max"
3335  << af_params.s[i] << af_params.D[i] << sw2 << sw2D << sqb << s_max;
3336 
3337  // Grid for moving adaptive FEM for faster sedimentation
3338 
3339  QVector< double > xbvec;
3340  xbvec.clear();
3341  xbvec.reserve( Nx );
3342  xbvec << af_params.current_meniscus;
3343 
3344  for ( int j = 1; j < Nx; j++ )
3345  {
3346  double dval = 0.1 * exp( sw2D * ( sq( 0.5 *
3347  ( xA[ j - 1 ] + xA[ j ] ) ) - sqb ) );
3348 
3349  double alpha = s_rat * ( 1.0 - dval ) + dval;
3350  xbvec << ( pow( xA[ j - 1 ], alpha ) *
3351  pow( xA[ j ], ( 1.0 - alpha ) ) );
3352 //if((j+1)==Nx) xbvec[j]=af_params.current_bottom;
3353 if(j<4 || j==(Nx/2) || (j+4)>Nx)
3354 DbgLv(2) << "RSA: smin>0:GlStf: j" << j << "dval alpha xbj xAj"
3355  << dval << alpha << xbvec[j] << xA[j];
3356  }
3357 
3358  double* xb = xbvec.data();
3359 
3360  GlobalStiff( xb, CA[ i ], CB[ i ], af_params.D[ i ], sw2 );
3361 DbgLv(2) << "RSA: smin>0:GlStf: CA[i]:" << CA[i][0][0] << CA[i][1][0]
3362  << CA[i][2][Nx-1];
3363 DbgLv(2) << "RSA: smin>0:GlStf: CB[i]:" << CB[i][0][0] << CB[i][1][0]
3364  << CB[i][2][Nx-1];
3365  }
3366  }
3367 
3368  else if ( s_max < 0) // all components floating
3369  {
3370  DbgErr() << "all components floating, not implemented yet";
3371  return -1;
3372  }
3373 
3374  else // sedimentation and floating mixed
3375  {
3376  DbgErr() << "sedimentation and floating mixed, suppose use "
3377  "fixed grid!";
3378  return -1;
3379  }
3380  } // moving mesh
3381  } // acceleration
3382 
3383  // Initial condition
3384  double** C0; // C[m][j]: current/next concentration of m-th component at x_j
3385  double** C1; // C[0...Ms-1][0....Nx-1]:
3386 
3387 DbgLv(2) << "RSA:_ra2:(6) Nx" << Nx << "x size" << x.size();
3388  US_AstfemMath::initialize_2d( Mcomp, Nx, &C0 );
3389  US_AstfemMath::initialize_2d( Mcomp, Nx, &C1 );
3390 
3391  // Here we need the interpolate the initial partial
3392  // concentration onto new grid x[j]
3393 
3394  for( int ii = 0; ii < Mcomp; ii++ )
3395  {
3396  // Interpolate the given C_init vector on the new C0 grid
3397  US_AstfemMath::interpolate_C0( C_init[ ii ], C0[ ii ], x );
3398  }
3399 
3400  // Total concentration at current and next time step
3401 DbgLv(2) << "RSA: newX3 Nx" << Nx << "C0[00] C0[m0] C0[0n] C0[mn]"
3402  << C0[0][0] << C0[Mcomp-1][0] << C0[0][Nx-1] << C0[Mcomp-1][Nx-1];
3403  QVector< double > CT0vec( Nx );
3404  QVector< double > CT1vec( Nx );
3405  QVector< double > rhVec ( Nx );
3406  double* CT0 = CT0vec.data();
3407  double* CT1 = CT1vec.data();
3408 
3409  for ( int jj = 0; jj < Nx; jj++ )
3410  {
3411  CT0[ jj ] = 0.0;
3412 
3413  for ( int ii = 0; ii < Mcomp; ii++ )
3414  CT0[ jj ] += C0[ ii ][ jj ];
3415 
3416  CT1[ jj ] = CT0[ jj ];
3417  }
3418 DbgLv(2) << "RSA: newX3 CT0 CTn" << CT1[0] << CT1[Nx-1];
3419 
3420  // Time evolution
3421  double* right_hand_side = rhVec.data();
3422 #ifndef NO_DB
3423  int stepinc = 1000;
3424  int stepmax = ( Nt + 2 ) / stepinc + 1;
3425  bool repprog = stepmax > 1;
3426 
3427  if ( repprog )
3428  {
3429  emit current_component( -stepmax );
3430  emit current_component( 0 );
3431  }
3432 #endif
3433 DbgLv(2) << "RSA: (5)AP.start_om2t" << af_params.start_om2t;
3434  last_time = af_params.start_time;
3435 
3436  for ( int kkk = 0; kkk < Nt + 2; kkk += 2 ) // two steps in together
3437  {
3438  double rpm_current = rpm_start +
3439  ( rpm_stop - rpm_start ) * ( kkk + 0.5 ) / Nt;
3440 
3441 #ifndef NO_DB
3442  emit current_speed( (int) rpm_current);
3443 #endif
3444 
3445  simscan.time = af_params.start_time + kkk * af_params.dt;
3446  simscan.rpm = (int) rpm_current;
3447  w2t_integral += ( ( simscan.time - last_time )
3448  * sq( rpm_current * M_PI / 30 ) );
3449  last_time = simscan.time;
3450  simscan.omega_s_t = w2t_integral;
3451 DbgLv(2) << "TMS:RSA:ra: time omegast" << simscan.time << simscan.omega_s_t
3452  << "step-scan" << simdata.scan.size();
3453 
3454  simscan.conc.clear();
3455  simscan.conc.reserve( Nx );
3456 DbgLv(2) << "RSA: Nx CT0size" << Nx << CT0vec.size() << "step-scan"
3457  << simdata.scan.size() << "rss-now" << US_Memory::rss_now();
3458 
3459  for ( int jj = 0; jj < Nx; jj++ )
3460  simscan.conc.append( CT0[ jj ] );
3461 DbgLv(2) << "TMS:RSA:ra: kkk" << kkk << "CT0[0] CT0[n]"
3462  << CT0[0] << CT0[Nx-1] << "accel fixedGrid" << accel << fixedGrid;
3463 
3464  simdata.scan .append( simscan );
3465 
3466  // First half step of sedimentation:
3467 
3468  if ( accel ) // need to reconstruct CA and CB by linear interpolation
3469  {
3470  double dval = sq( rpm_current / rpm_stop );
3471 
3472  for ( int i = 0; i < Mcomp; i++ )
3473  {
3474  for ( int j1 = 0; j1 < 3; j1++ )
3475  {
3476  for ( int j2 = 0; j2 < Nx; j2++ )
3477  {
3478  CA[ i ][ j1 ][ j2 ] = CA1[ i ][ j1 ][ j2 ] +
3479  dval * ( CA2[ i ][ j1 ][ j2 ] -
3480  CA1[ i ][ j1 ][ j2 ] );
3481 
3482  CB[ i ][ j1 ][ j2 ] = CB1[ i ][ j1 ][ j2 ] +
3483  dval * ( CB2[ i ][ j1 ][ j2 ] -
3484  CB1[ i ][ j1 ][ j2 ] );
3485  }
3486  }
3487  }
3488  }
3489 
3490  if ( accel || fixedGrid ) // For fixed grid
3491  {
3492  for ( int i = 0; i < Mcomp; i++ )
3493  {
3494  right_hand_side[ 0 ] = - CB[ i ][ 1 ][ 0 ] * C0[ i ][ 0 ]
3495  - CB[ i ][ 2 ][ 0 ] * C0[ i ][ 1 ];
3496 
3497  for ( int j = 1; j < Nx - 1; j++ )
3498  {
3499  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 1 ]
3500  - CB[ i ][ 1 ][ j ] * C0[ i ][ j ]
3501  - CB[ i ][ 2 ][ j ] * C0[ i ][ j + 1 ];
3502  }
3503 
3504  int j = Nx - 1;
3505 
3506  right_hand_side[ j ] = -CB[ i ][ 0 ][ j ] * C0[ i ][ j - 1 ]
3507  - CB[ i ][ 1 ][ j ] * C0[ i ][ j ];
3508 
3509  US_AstfemMath::tridiag( CA[ i ][ 0 ], CA[ i ][ 1 ], CA[ i ][ 2 ],
3510  right_hand_side, C1[ i ], Nx );
3511  }
3512  }
3513 
3514  else // Moving grid
3515  {
3516 //DbgLv(2) << "TMS:RSA:ra: (0)MovGrid: C0[00] C0[0n] C0[m0] C0[1n]"
3517 // << C0[0][0] << C0[0][Nx-1] << C0[Mcomp-1][0] << C0[Mcomp-1][Nx-1];
3518 //DbgLv(2) << "TMS:RSA:ra: (0)MovGrid: C1[00] C1[0n] C1[m0] C1[1n]"
3519 // << C1[0][0] << C1[0][Nx-1] << C1[Mcomp-1][0] << C1[Mcomp-1][Nx-1];
3520  for ( int i = 0; i < Mcomp; i++ )
3521  {
3522  // Calculate the right hand side vector
3523  right_hand_side[ 0 ] = - CB[ i ][ 2 ][ 0 ] * C0[ i ][ 0 ]
3524  - CB[ i ][ 3 ][ 0 ] * C0[ i ][ 1 ];
3525 
3526  right_hand_side[ 1 ] = - CB[ i ][ 1 ][ 1 ] * C0[ i ][ 0 ]
3527  - CB[ i ][ 2 ][ 1 ] * C0[ i ][ 1 ]
3528  - CB[ i ][ 3 ][ 1 ] * C0[ i ][ 2 ];
3529 
3530  for ( int j = 2; j < Nx - 1; j++ )
3531  {
3532  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 2 ]
3533  - CB[ i ][ 1 ][ j ] * C0[ i ][ j - 1 ]
3534  - CB[ i ][ 2 ][ j ] * C0[ i ][ j ]
3535  - CB[ i ][ 3 ][ j ] * C0[ i ][ j + 1 ];
3536  }
3537 
3538  int j = Nx - 1;
3539  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 2 ]
3540  - CB[ i ][ 1 ][ j ] * C0[ i ][ j - 1 ]
3541  - CB[ i ][ 2 ][ j ] * C0[ i ][ j ];
3542 
3543 DbgLv(2) << "TMS:RSA:ra: MovGrid: i" << i << "CB[0]" << CB[i][2][0]
3544  << CB[i][3][0] << CB[i][1][1] << CB[i][2][1] << CB[i][3][1];
3545 DbgLv(2) << "TMS:RSA:ra: MovGrid: i" << i << "CA[0]" << CA[i][0][0]
3546  << CA[i][1][0] << CA[i][0][0] << CA[i][1][0] << CA[i][2][0];
3547  US_AstfemMath::QuadSolver( CA[ i ][ 0 ], CA[ i ][ 1 ], CA[ i ][ 2 ],
3548  CA[ i ][ 3 ], right_hand_side, C1[ i ],
3549  Nx );
3550 DbgLv(2) << "TMS:RSA:ra: MovGrid: i" << i << "C1[i]" << C1[i][0]
3551  << C1[i][Nx-4] << C1[i][Nx-3] << C1[i][Nx-2] << C1[i][Nx-1];
3552 DbgLv(2) << "TMS:RSA:ra: MovGrid: i" << i << "C0[i]" << C0[i][0]
3553  << C0[i][Nx-4] << C0[i][Nx-3] << C0[i][Nx-2] << C0[i][Nx-1];
3554  }
3555 //DbgLv(2) << "TMS:RSA:ra: MovGrid: CA[000] rs[0]" << CA[0][0][0]
3556 // << right_hand_side[0] << "CB[000]" << CB[0][0][0];
3557 //DbgLv(2) << "TMS:RSA:ra: MovGrid: C1[00] C1[0n] C1[m0] C1[1n]"
3558 // << C1[0][0] << C1[0][Nx-1] << C1[Mcomp-1][0] << C1[Mcomp-1][Nx-1];
3559  }
3560 
3561  // Reaction part: instantaneous reaction at each node
3562  //
3563  // instantaneous reaction at each node
3564  // [C1]=ReactionOneStep_inst(C1);
3565  //
3566  // Finite reaction rate: linear interpolation of instantaneous reaction
3567 
3568 DbgLv(2) << "TMS:RSA:ra: C1[10] C1[11] C1[1k] C1[1n]"
3569  << C1[1][0] << C1[1][1] << C1[1][Nx-2] << C1[1][Nx-1];
3570  ReactionOneStep_Euler_imp( Nx, C1, 2 * af_params.dt );
3571 DbgLv(2) << "TMS:RSA:ra(2): C1[10] C1[11] C1[1k] C1[1n]"
3572  << C1[1][0] << C1[1][1] << C1[1][Nx-2] << C1[1][Nx-1];
3573 
3574  // For next half time-step in SNI operator splitting scheme
3575 
3576  for ( int j = 0; j < Nx; j++ )
3577  {
3578  CT1[ j ] = 0.0;
3579 
3580  for ( int i = 0; i < Mcomp; i++ )
3581  {
3582  CT1[ j ] += C1[ i ][ j ];
3583  C0[ i][ j ] = C1[ i ][ j ];
3584  }
3585 
3586  CT0[ j ] = CT1[ j ];
3587  }
3588 
3589  // 2nd half step of sedimentation:
3590 
3591  rpm_current = rpm_start + ( rpm_stop - rpm_start ) * ( kkk + 1.5 ) / Nt;
3592 
3593  if ( accel ) // Need to reconstruct CA and CB by linear interpolation
3594  {
3595  double dval = sq( rpm_current / rpm_stop );
3596 
3597  for ( int i = 0; i < Mcomp; i++ )
3598  {
3599  for ( int j1 = 0; j1 < 3; j1++ )
3600  {
3601  for ( int j2 = 0; j2 < Nx; j2++ )
3602  {
3603  CA[ i][ j1 ][ j2 ] = CA1[ i ][ j1 ][ j2 ] +
3604  dval * ( CA2[ i ][ j1 ][ j2 ] -
3605  CA1[ i ][ j1 ][ j2 ] );
3606 
3607  CB[ i][ j1 ][ j2 ] = CB1[ i ][ j1 ][ j2 ] +
3608  dval * ( CB2[ i ][ j1 ][ j2 ] -
3609  CB1[ i ][ j1 ][ j2 ] ) ;
3610  }
3611  }
3612  }
3613  }
3614 
3615  if ( accel || fixedGrid ) // For fixed grid
3616  {
3617  for ( int i = 0; i < Mcomp; i++ )
3618  {
3619  right_hand_side[ 0 ] = - CB[ i ][ 1 ][ 0 ] * C0[ i ][ 0 ]
3620  - CB[ i ][ 2 ][ 0 ] * C0[ i ][ 1 ];
3621 
3622  for ( int j = 1; j < Nx - 1; j++ )
3623  {
3624  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 1 ]
3625  - CB[ i ][ 1 ][ j ] * C0[ i ][ j ]
3626  - CB[ i ][ 2 ][ j ] * C0[ i ][ j + 1 ];
3627  }
3628 
3629  int j = Nx - 1;
3630  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 1 ]
3631  - CB[ i ][ 1 ][ j ] * C0[ i ][ j ];
3632 
3633  US_AstfemMath::tridiag( CA[ i ][ 0 ], CA[ i ][ 1 ], CA[ i ][ 2 ],
3634  right_hand_side, C1[ i ], Nx );
3635  }
3636  }
3637  else // Moving grid
3638  {
3639  for ( int i = 0; i < Mcomp; i++ )
3640  {
3641  // Calculate the right hand side vector //
3642  right_hand_side[ 0 ] = - CB[ i ][ 2 ][ 0 ] * C0[ i ][ 0 ]
3643  - CB[ i ][ 3 ][ 0 ] * C0[ i ][ 1 ];
3644 
3645  right_hand_side[ 1 ] = - CB[ i ][ 1 ][ 1 ] * C0[ i ][ 0 ]
3646  - CB[ i ][ 2 ][ 1 ] * C0[ i ][ 1 ]
3647  - CB[ i ][ 3 ][ 1 ] * C0[ i ][ 2 ];
3648 
3649  for ( int j = 2; j < Nx - 1; j++ )
3650  {
3651  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 2 ]
3652  - CB[ i ][ 1 ][ j ] * C0[ i ][ j - 1 ]
3653  - CB[ i ][ 2 ][ j ] * C0[ i ][ j ]
3654  - CB[ i ][ 3 ][ j ] * C0[ i ][ j + 1 ];
3655  }
3656 
3657  int j = Nx - 1;
3658  right_hand_side[ j ] = - CB[ i ][ 0 ][ j ] * C0[ i ][ j - 2 ]
3659  - CB[ i ][ 1 ][ j ] * C0[ i ][ j - 1 ]
3660  - CB[ i ][ 2 ][ j ] * C0[ i ][ j ];
3661 
3662  US_AstfemMath::QuadSolver( CA[ i ][ 0 ], CA[ i ][ 1 ],
3663  CA[ i ][ 2 ], CA[ i ][ 3 ],
3664  right_hand_side, C1[ i ], Nx );
3665  }
3666 //DbgLv(2) << "TMS:RSA:ra: MovGrid: C1[00] C1[0n] C1[m0] C1[1n]"
3667 // << C1[0][0] << C1[0][Nx-1] << C1[Mcomp-1][0] << C1[Mcomp-1][Nx-1];
3668  }
3669 
3670  // End of 2nd half step of sedimentation
3671 
3672  // For next 2 time steps
3673 
3674  for ( int j = 0; j < Nx; j++ )
3675  {
3676  CT1[ j ] = 0.0;
3677 
3678  for ( int i = 0; i < Mcomp; i++ )
3679  {
3680  CT1[ j ] += C1[ i ][ j ];
3681  C0[ i ][ j ] = C1[ i ][ j ];
3682  }
3683 
3684  CT0[ j ] = CT1[ j ];
3685  }
3686 
3687 #ifndef NO_DB
3688  //if ( show_movie && (kkk%8) == 0 )
3689  if ( show_movie )
3690  {
3691  if ( stopFlag ) break;
3692 
3693  emit new_scan( &x, CT0 );
3694  emit new_time( simscan.time );
3695  qApp->processEvents();
3696  }
3697 
3698  if ( repprog && ( ( kkk + 1 ) & stepinc ) == 0 )
3699  {
3700  emit current_component( ( kkk + 1 ) / stepinc );
3701  }
3702 #endif
3703 
3704  } // time loop
3705 
3706 #ifndef NO_DB
3707  emit new_time( simscan.time );
3708  qApp->processEvents();
3709 #endif
3710 
3711  for ( int i = 0; i < Mcomp; i++ )
3712  {
3713  C_init[ i ].radius .clear();
3714  C_init[ i ].concentration.clear();
3715  C_init[ i ].radius .reserve( Nx );
3716  C_init[ i ].concentration.reserve( Nx );
3717 
3718  for ( int j = 0; j < Nx; j++ )
3719  {
3720  C_init[ i ].radius .append( xA[ j ] );
3721  C_init[ i ].concentration .append( C1[ i ][ j ] );
3722  }
3723  }
3724 
3725  US_AstfemMath::clear_2d(Mcomp, C0);
3726  US_AstfemMath::clear_2d(Mcomp, C1);
3727  US_AstfemMath::clear_3d(Mcomp, 4, CA);
3728  US_AstfemMath::clear_3d(Mcomp, 4, CB);
3729 
3730  if ( accel ) // then we have acceleration
3731  {
3732  US_AstfemMath::clear_3d( Mcomp, 3, CA1 );
3733  US_AstfemMath::clear_3d( Mcomp, 3, CB1 );
3734  US_AstfemMath::clear_3d( Mcomp, 3, CA2 );
3735  US_AstfemMath::clear_3d( Mcomp, 3, CB2 );
3736  }
3737 
3738  return 0;
3739 }
3740 
3741 void US_Astfem_RSA::GlobalStiff( double* xb, double** ca,
3742  double** cb, double D, double sw2 )
3743 {
3744  // 4: global stifness matrix
3745 
3746  // function [CA, CB]=4(x, xb, dt, D, sw2)
3747 
3748  double*** Stif = NULL;
3749  double vx[ 8 ];
3750 
3751  US_AstfemMath::initialize_3d( Nx, 6, 2, &Stif );
3752 
3753  // 1st elem
3754  vx[ 0 ] = x [ 0 ];
3755  vx[ 1 ] = x [ 1 ];
3756  vx[ 2 ] = x [ 0 ];
3757  vx[ 3 ] = x [ 1 ];
3758  vx[ 4 ] = x [ 2 ];
3759  vx[ 5 ] = xb[ 1 ];
3760  US_AstfemMath::IntQT1( vx, D, sw2, Stif[ 0 ], af_params.dt );
3761 
3762  // elems in middle
3763  for ( int i = 1; i < Nx - 2; i++ )
3764  {
3765  vx[ 0 ] = x [ i - 1 ];
3766  vx[ 1 ] = x [ i ];
3767  vx[ 2 ] = x [ i + 1 ];
3768  vx[ 3 ] = x [ i ];
3769  vx[ 4 ] = x [ i + 1 ];
3770  vx[ 5 ] = x [ i + 2 ];
3771  vx[ 6 ] = xb[ i ];
3772  vx[ 7 ] = xb[ i + 1 ];
3773  US_AstfemMath::IntQTm( vx, D, sw2, Stif[ i ], af_params.dt );
3774  }
3775 
3776  // 2nd last elems
3777  vx[ 0 ] = x [ Nx - 3 ];
3778  vx[ 1 ] = x [ Nx - 2 ];
3779  vx[ 2 ] = x [ Nx - 1 ];
3780  vx[ 3 ] = x [ Nx - 2 ];
3781  vx[ 4 ] = x [ Nx - 1 ];
3782  vx[ 5 ] = xb[ Nx - 2 ];
3783  vx[ 6 ] = xb[ Nx - 1 ];
3784 
3785  US_AstfemMath::IntQTn2( vx, D, sw2, Stif[ Nx - 2 ], af_params.dt );
3786 
3787  // last elems
3788  vx[ 0 ] = x [ Nx - 2 ];
3789  vx[ 1 ] = x [ Nx - 1 ];
3790  vx[ 2 ] = x [ Nx - 1 ];
3791  vx[ 3 ] = xb[ Nx - 1 ];
3792  US_AstfemMath::IntQTn1 ( vx, D, sw2, Stif[ Nx - 1 ], af_params.dt );
3793 
3794  // assembly into global stiffness matrix
3795  ca[ 0 ][ 0 ] = 0.0;
3796  ca[ 1 ][ 0 ] = Stif[ 0 ][ 2 ][ 0 ];
3797  ca[ 2 ][ 0 ] = Stif[ 0 ][ 3 ][ 0 ];
3798  ca[ 3 ][ 0 ] = Stif[ 0 ][ 4 ][ 0 ];
3799 
3800  cb[ 0 ][ 0 ] = 0.0;
3801  cb[ 1 ][ 0 ] = 0.0;
3802  cb[ 2 ][ 0 ] = Stif[ 0 ][ 0 ][ 0 ];
3803  cb[ 3 ][ 0 ] = Stif[ 0 ][ 1 ][ 0 ];
3804 
3805  // i=2
3806  ca[ 0 ][ 1 ] = Stif[ 0 ][ 2 ][ 1 ];
3807  ca[ 1 ][ 1 ] = Stif[ 0 ][ 3 ][ 1 ] + Stif[ 1 ][ 3 ][ 0 ];
3808  ca[ 2 ][ 1 ] = Stif[ 0 ][ 4 ][ 1 ] + Stif[ 1 ][ 4 ][ 0 ];
3809  ca[ 3 ][ 1 ] = Stif[ 1 ][ 5 ][ 0 ];
3810 
3811  cb[ 0 ][ 1 ] = 0.0;
3812  cb[ 1 ][ 1 ] = Stif[ 0 ][ 0 ][ 1 ] + Stif[ 1 ][ 0 ][ 0 ];
3813  cb[ 2 ][ 1 ] = Stif[ 0 ][ 1 ][ 1 ] + Stif[ 1 ][ 1 ][ 0 ];
3814  cb[ 3 ][ 1 ] = Stif[ 1 ][ 2 ][ 0 ];
3815 
3816  // i: middle
3817  for ( int i = 2; i < Nx - 2; i++ )
3818  {
3819  ca[ 0 ][ i ] = Stif[ i - 1 ][ 3 ][ 1 ];
3820  ca[ 1 ][ i ] = Stif[ i - 1 ][ 4 ][ 1 ] + Stif[ i ][ 3 ][ 0 ];
3821  ca[ 2 ][ i ] = Stif[ i - 1 ][ 5 ][ 1 ] + Stif[ i ][ 4 ][ 0 ];
3822  ca[ 3 ][ i ] = Stif[ i ][ 5 ][ 0 ];
3823 
3824  cb[ 0 ][ i ] = Stif[ i - 1 ][ 0 ][ 1 ];
3825  cb[ 1 ][ i ] = Stif[ i - 1 ][ 1 ][ 1 ] + Stif[ i ][ 0 ][ 0 ];
3826  cb[ 2 ][ i ] = Stif[ i - 1 ][ 2 ][ 1 ] + Stif[ i ][ 1 ][ 0 ];
3827  cb[ 3 ][ i ] = Stif[ i ][ 2 ][ 0 ];
3828  }
3829 
3830  // i=n
3831  int i = Nx - 2;
3832  ca[ 0 ][ i ] = Stif[ i - 1 ][ 3 ][ 1 ];
3833  ca[ 1 ][ i ] = Stif[ i - 1 ][ 4 ][ 1 ] + Stif[ i ][ 3 ][ 0 ];
3834  ca[ 2 ][ i ] = Stif[ i - 1 ][ 5 ][ 1 ] + Stif[ i ][ 4 ][ 0 ];
3835  ca[ 3 ][ i ] = 0.0;
3836 
3837  cb[ 0 ][ i ] = Stif[ i - 1 ][ 0 ][ 1 ];
3838  cb[ 1 ][ i ] = Stif[ i - 1 ][ 1 ][ 1 ] + Stif[ i ][ 0 ][ 0 ];
3839  cb[ 2 ][ i ] = Stif[ i - 1 ][ 2 ][ 1 ] + Stif[ i ][ 1 ][ 0 ];
3840  cb[ 3 ][ i ] = Stif[ i ][ 2 ][ 0 ];
3841 
3842  // i=n+1
3843  i = Nx - 1;
3844  ca[ 0 ][ i ] = Stif[ i - 1 ][ 3 ][ 1 ];
3845  ca[ 1 ][ i ] = Stif[ i - 1 ][ 4 ][ 1 ] + Stif[ i ][ 2 ][ 0 ];
3846  ca[ 2 ][ i ] = 0.0;
3847  ca[ 3 ][ i ] = 0.0;
3848 
3849  cb[ 0 ][ i ] = Stif[ i - 1 ][ 0 ][ 1 ];
3850  cb[ 1 ][ i ] = Stif[ i - 1 ][ 1 ][ 1 ] + Stif[ i ][ 0 ][ 0 ];
3851  cb[ 2 ][ i ] = Stif[ i - 1 ][ 2 ][ 1 ] + Stif[ i ][ 1 ][ 0 ];
3852  cb[ 3 ][ i ] = 0.0;
3853 
3854  US_AstfemMath::clear_3d( Nx, 6, Stif );
3855 }
3856 
3857 void US_Astfem_RSA::load_mfem_data( US_DataIO::RawData& edata,
3858  US_AstfemMath::MfemData& fdata )
3859 {
3860  int nscan = edata.scanCount();
3861  int nconc = edata.pointCount();
3862 DbgLv(2) << "RSA:f nscan nconc" << nscan << nconc;
3863 
3864  fdata.id = edata.description;
3865  fdata.cell = edata.cell;
3866  fdata.scan .resize( nscan );
3867 // fdata.radius.resize( nconc );
3868  fdata.radius = edata.xvalues;
3869 DbgLv(2) << "RSA:f r0 rn" << fdata.radius[0] << fdata.radius[nconc-1];
3870 
3871  for ( int ii = 0; ii < nscan; ii++ )
3872  {
3873  US_AstfemMath::MfemScan* fscan = &fdata.scan[ ii ];
3874 
3875  fscan->temperature = edata.scanData[ ii ].temperature;
3876  fscan->rpm = edata.scanData[ ii ].rpm;
3877  fscan->time = edata.scanData[ ii ].seconds;
3878  fscan->omega_s_t = edata.scanData[ ii ].omega2t;
3879  fscan->conc = edata.scanData[ ii ].rvalues;
3880 if(ii<3)
3881 DbgLv(2) << "RSA:f ii c0 cn" << ii << fscan->conc[0] << fscan->conc[nconc-1]
3882  << "d0 dn" << edata.scanData[ii].rvalues[0] << edata.scanData[ii].rvalues[nconc-1];
3883 #if 1
3884  }
3885 #endif
3886 #if 0
3887  fscan->conc.resize( nconc );
3888 
3889  for ( int jj = 0; jj < nconc; jj++ )
3890  {
3891  fscan->conc[ jj ] = edata.value( ii, jj );
3892  }
3893  }
3894 
3895  for ( int jj = 0; jj < nconc; jj++ )
3896  {
3897  fdata.radius[ jj ] = edata.radius( jj );
3898  }
3899 #endif
3900 DbgLv(2) << "RSA:f sc0 temp" << fdata.scan[0].temperature;
3901 DbgLv(2) << "RSA:e sc0 temp" << edata.scanData[0].temperature;
3902 }
3903 
3904 void US_Astfem_RSA::store_mfem_data( US_DataIO::RawData& edata,
3905  US_AstfemMath::MfemData& fdata )
3906 {
3907  int nscan = fdata.scan.size();
3908  int nconc = fdata.radius.size();
3909  int escan = edata.scanCount();
3910  int econc = edata.pointCount();
3911 DbgLv(2) << "RSA:st_md: nscan nconc" << nscan << nconc;
3912 DbgLv(2) << "RSA:st_md: escan econc" << escan << econc;
3913 
3914  if ( escan != nscan )
3915  edata.scanData.resize( nscan );
3916 
3917  edata.description = fdata.id;
3918  edata.cell = fdata.cell;
3919  edata.xvalues = fdata.radius;
3920 
3921  for ( int ii = 0; ii < nscan; ii++ )
3922  {
3923  US_AstfemMath::MfemScan* fscan = &fdata.scan [ ii ];
3924  US_DataIO::Scan* escan = &edata.scanData[ ii ];
3925 
3926  escan->temperature = fscan->temperature;
3927  escan->rpm = fscan->rpm;
3928  escan->seconds = fscan->time;
3929  escan->omega2t = fscan->omega_s_t;
3930  escan->plateau = fdata.radius[ nconc - 1 ];
3931  escan->rvalues = fscan->conc;
3932  }
3933 
3934 DbgLv(2) << "RSA:o-f sc0 temp" << fdata.scan[0].temperature;
3935 DbgLv(2) << "RSA:o-e sc0 temp" << edata.scanData[0].temperature;
3936 }
Generated on Sun Jun 21 2015 11:40:41 for UltraScan III by   doxygen 1.8.3.1-20130324