UltraScan III
2dsa_worker.cpp
Go to the documentation of this file.
1 #include "us_mpi_analysis.h"
2 #include "us_constants.h"
3 #include "us_math2.h"
4 #include "us_astfem_rsa.h"
5 #include "us_simparms.h"
6 
7 void US_MPI_Analysis::_2dsa_worker( void )
8 {
9  bool repeat_loop = true;
10  MPI_Job job;
11  MPI_Status status;
12 
13  // Use 4 here because the master will be reading 4 with the
14  // same instruction when reading ::READY or ::RESULTS.
15  int x[ 4 ];
16 
17  while ( repeat_loop )
18  {
19  MPI_Send( x, // Basically don't care
20  4,
21  MPI_INT,
22  MPI_Job::MASTER,
23  MPI_Job::READY,
24  my_communicator ); // let master know we are ready
25 //if(my_rank==1)
26 DbgLv(1) << "w:" << my_rank << ": ready sent";
27 
28  // Blocking -- Wait for instructions
29  MPI_Recv( &job, // get masters' response
30  sizeof( job ),
31  MPI_BYTE,
32  MPI_Job::MASTER,
33  MPI_Job::TAG0,
34  my_communicator,
35  &status ); // status not used
36 //if(my_rank==1)
37 DbgLv(1) << "w:" << my_rank << ": job_recvd length" << job.length
38  << "command" << job.command;
39 
40  meniscus_value = job.meniscus_value;
41  int offset = job.dataset_offset;
42  int dataset_count = job.dataset_count;
43  int job_length = job.length;
44 DbgLv(1) << "w:" << my_rank << ": offs cnt" << offset << dataset_count;
45 
46  data_sets[ offset ]->run_data.meniscus = meniscus_value;
47  data_sets[ offset ]->simparams.meniscus = meniscus_value;
48 
49  switch( job.command )
50  {
51 
52  case MPI_Job::PROCESS: // Process solutes
53  {
54  US_SolveSim::Simulation simulation_values;
55  simulation_values.noisflag =
56  parameters[ "tinoise_option" ].toInt() > 0 ? 1 : 0;
57  simulation_values.noisflag +=
58  parameters[ "rinoise_option" ].toInt() > 0 ? 2 : 0;
59  simulation_values.dbg_level = dbg_level;
60  simulation_values.dbg_timing = dbg_timing;
61 
62 //DbgLv(1) << "w:" << my_rank << ": sols size" << job.length;
63 //if(my_rank==1)
64 DbgLv(1) << "w:" << my_rank << ": sols size" << job.length;
65  simulation_values.solutes.resize( job.length );
66 
67  MPI_Recv( simulation_values.solutes.data(), // Get solutes
68  job_length * solute_doubles,
69  MPI_DOUBLE,
70  MPI_Job::MASTER,
71  MPI_Job::TAG0,
72  my_communicator,
73  &status );
74 
75  max_rss();
76 //*DEBUG*
77 //if(dbg_level>0 && my_rank==1)
78 //if(my_rank==1)
79 {
80  int nn = simulation_values.solutes.size() - 1;
81  int mm = nn/2;
82  DbgLv(1) << "w:" << my_rank << ": offs dscnt" << offset << dataset_count
83  << "vbar s20wc bott"
84  << data_sets[offset]->vbar20
85  << data_sets[offset]->s20w_correction
86  << data_sets[offset]->centerpiece_bottom;
87  DbgLv(1) << "w:" << my_rank << ": sol0 solm soln"
88  << simulation_values.solutes[0].s << simulation_values.solutes[0].k
89  << simulation_values.solutes[mm].s << simulation_values.solutes[mm].k
90  << simulation_values.solutes[nn].s << simulation_values.solutes[nn].k;
91 }
92 //*DEBUG*
93 
94  calc_residuals( offset, dataset_count, simulation_values );
95 
96  // Tell master we are sending back results
97  int size[ 4 ] = { simulation_values.solutes.size(),
98  simulation_values.ti_noise.size(),
99  simulation_values.ri_noise.size(),
100  max_rss() };
101 
102 DbgLv(1) << "w:" << my_rank << ": result sols size" << size[0]
103  << "max_rss" << size[ 3 ];
104 //*DEBUG*
105 if(dbg_level==0 && my_rank==1) {
106 DbgLv(1) << "w:" << my_rank << ": result sols size" << size[0]
107  << "nsscan" << simulation_values.sim_data.scanCount();
108 }
109 //*DEBUG*
110  MPI_Send( size,
111  4,
112  MPI_INT,
113  MPI_Job::MASTER,
114  MPI_Job::RESULTS,
115  my_communicator );
116 
117  // Send back to master all of simulation_values
118  MPI_Send( simulation_values.solutes.data(),
119  simulation_values.solutes.size() * solute_doubles,
120  MPI_DOUBLE,
121  MPI_Job::MASTER,
122  MPI_Job::TAG0,
123  my_communicator );
124 
125  MPI_Send( &simulation_values.variance,
126  1,
127  MPI_DOUBLE,
128  MPI_Job::MASTER,
129  MPI_Job::TAG0,
130  my_communicator );
131 
132  MPI_Send( simulation_values.variances.data(),
133  dataset_count,
134  MPI_DOUBLE,
135  MPI_Job::MASTER,
136  MPI_Job::TAG0,
137  my_communicator );
138 
139  MPI_Send( simulation_values.ti_noise.data(),
140  simulation_values.ti_noise.size(),
141  MPI_DOUBLE,
142  MPI_Job::MASTER,
143  MPI_Job::TAG0,
144  my_communicator );
145 
146  MPI_Send( simulation_values.ri_noise.data(),
147  simulation_values.ri_noise.size(),
148  MPI_DOUBLE,
149  MPI_Job::MASTER,
150  MPI_Job::TAG0,
151  my_communicator );
152  }
153 
154  break;
155 
156  case MPI_Job::NEWDATA: // Reset data for Monte Carlo or global fit
157  {
158  int mc_iter = job.solution;
159 
160  //if ( dataset_count > 0 && mc_iter < 4 )
161  if ( is_global_fit && mc_iter < 3 && my_rank < 3 )
162  { // For global fits, check the memory requirements
163  long memused = max_rss();
164  long memdata = job_length * sizeof( double );
165  int grid_reps = qMax( parameters[ "uniform_grid" ].toInt(), 1 );
166  double s_pts = 60.0;
167  double ff0_pts = 60.0;
168  if ( parameters.contains( "s_grid_points" ) )
169  s_pts = parameters[ "s_grid_points" ].toDouble();
170  else if ( parameters.contains( "s_resolution" ) )
171  s_pts = parameters[ "s_resolution" ].toDouble() * grid_reps;
172  if ( parameters.contains( "ff0_grid_points" ) )
173  ff0_pts = parameters[ "ff0_grid_points" ].toDouble();
174  else if ( parameters.contains( "ff0_resolution" ) )
175  ff0_pts = parameters[ "ff0_resolution" ].toDouble() * grid_reps;
176  int nsstep = (int)( s_pts );
177  int nkstep = (int)( ff0_pts );
178  int maxsols = nsstep * nkstep;
179  long memamatr = memdata * maxsols;
180  long membmatr = memdata;
181  long memneed = memdata + memamatr + membmatr;
182  const double mb_bytes = ( 1024. * 1024. );
183  const double gb_bytes = ( mb_bytes * 1024. );
184  double gb_need = (double)memneed / gb_bytes;
185  gb_need = qRound( gb_need * 1000.0 ) * 0.001;
186  double gb_used = (double)memused / mb_bytes;
187  gb_used = qRound( gb_used * 1000.0 ) * 0.001;
188  long pgavail = sysconf( _SC_PHYS_PAGES );
189  long pgsize = sysconf( _SC_PAGE_SIZE );
190  long memavail = pgavail * pgsize;
191  double gb_avail = (double)memavail / gb_bytes;
192  gb_avail = qRound( gb_avail * 1000.0 ) * 0.001;
193  long pgcurav = sysconf( _SC_AVPHYS_PAGES );
194  long memcurav = pgcurav * pgsize;
195  double gb_curav = (double)memcurav / gb_bytes;
196  gb_curav = qRound( gb_curav * 1000.0 ) * 0.001;
197 
198  qDebug() << "++ Worker" << my_rank << ": MC iteration"
199  << mc_iter << ": Memory Profile :"
200  << "\n Maximum memory used to this point" << memused
201  << "\n Composite data memory needed" << memdata
202  << "\n Maximum subgrid solute count" << maxsols
203  << "\n NNLS A matrix memory needed" << memamatr
204  << "\n NNLS B matrix memory needed" << membmatr
205  << "\n Total memory (GB) used" << gb_used
206  << "\n Total memory (GB) needed" << gb_need
207  << "\n Total memory (GB) available" << gb_avail
208  << "\n Memory (GB) currently available" << gb_curav;
209  }
210 
211  mc_data.resize( job_length );
212 
213  if ( mc_data.size() != job_length )
214  {
215  DbgLv(0) << "*ERROR* mc_data.size() job_length"
216  << mc_data.size() << job_length;
217  }
218 
219  MPI_Barrier( my_communicator );
220 
221 if(my_rank==1 || my_rank==11)
222 DbgLv(1) << "newD:" << my_rank << " scld/newdat rcv : offs dsknt"
223  << offset << dataset_count << "joblen" << job_length;
224 double dsum=0.0;
225  // This is a receive
226  MPI_Bcast( mc_data.data(),
227  job_length,
228  MPI_DOUBLE,
229  MPI_Job::MASTER,
230  my_communicator );
231 
232 
233  if ( is_global_fit && dataset_count == 1 )
234  { // For global update to scaled data, extra value is new ODlimit
235  job_length--;
236  data_sets[ offset ]->run_data.ODlimit = mc_data[ job_length ];
237 if( (my_rank==1||my_rank==11) )
238 DbgLv(1) << "newD:" << my_rank << ":offset ODlimit" << offset
239  << data_sets[ offset ]->run_data.ODlimit;
240  }
241 
242  int index = 0;
243 
244  for ( int ee = offset; ee < offset + dataset_count; ee++ )
245  {
246  US_DataIO::EditedData* edata = &data_sets[ ee ]->run_data;
247 
248  int scan_count = edata->scanCount();
249  int radius_points = edata->pointCount();
250 
251 int indxh=((scan_count/2)*radius_points)+(radius_points/2);
252  for ( int ss = 0; ss < scan_count; ss++ )
253  {
254  for ( int rr = 0; rr < radius_points; rr++, index++ )
255  {
256  edata->setValue( ss, rr, mc_data[ index ] );
257 dsum+=edata->value(ss,rr);
258 if( (my_rank==1||my_rank==11)
259  && (index<5 || index>(job_length-6) || (index>(indxh-4)&&index<(indxh+3))) )
260 DbgLv(1) << "newD:" << my_rank << ":index" << index << "edat" << edata->value(ss,rr)
261  << "ee" << ee;
262  }
263  }
264  }
265 if(my_rank==1 || my_rank==11)
266 DbgLv(1) << "newD:" << my_rank << " length index" << job_length << index
267  << "dsum" << dsum;
268  }
269 
270  break;
271 
272  default:
273  repeat_loop = false;
274  break;
275  } // switch
276  } // repeat_loop
277 }
278 
279 void US_MPI_Analysis::calc_residuals( int offset,
280  int dataset_count,
281  US_SolveSim::Simulation& simu_values )
282 {
283  count_calc_residuals++;
284 
285  US_SolveSim solvesim( data_sets, my_rank, false );
286 
287 //*DEBUG*
288 int dbglvsv=simu_values.dbg_level;
289 simu_values.dbg_level=(dbglvsv>1||my_rank==1)?dbglvsv:0;
290 //simu_values.dbg_level=(dbglvsv>0)?dbglvsv:0;
291 int nsoli=simu_values.solutes.size();
292 QVector< US_Solute > isols = simu_values.solutes;
293 //*DEBUG*
294 
295  solvesim.calc_residuals( offset, dataset_count, simu_values );
296 
297 //*DEBUG*
298 simu_values.dbg_level=dbglvsv;
299 if ( dbg_level > 0 && ( group_rank == 1 || group_rank == 11 ) )
300 //if ( group_rank == 1 || group_rank == 11 )
301 {
302  int nsolo=simu_values.solutes.size();
303  US_DataIO::EditedData* edat = &data_sets[offset]->run_data;
304  US_SolveSim::DataSet* dset = data_sets[offset];
305  US_DataIO::RawData* sdat = &simu_values.sim_data;
306  int nsc = edat->scanCount();
307  int nrp = edat->pointCount();
308  double d0 = edat->scanData[0].rvalues[0];
309  double d1 = edat->scanData[0].rvalues[1];
310  double dh = edat->scanData[nsc/2].rvalues[nrp/2];
311  double dm = edat->scanData[nsc-1].rvalues[nrp-2];
312  double dn = edat->scanData[nsc-1].rvalues[nrp-1];
313  DbgLv(1) << "w:" << my_rank << ":d(01hmn)" << d0 << d1 << dh << dm << dn;
314  double dt = 0.0;
315  for ( int ss=0;ss<nsc;ss++ )
316  for ( int rr=0;rr<nrp;rr++ ) dt += edat->scanData[ss].rvalues[rr];
317  DbgLv(1) << "w:" << my_rank << ":dtot" << dt;
318  double s0 = sdat->value(0,0);
319  double s1 = sdat->value(0,1);
320  double sh = sdat->value(nsc/2,nrp/2);
321  double sm = sdat->value(nsc-1,nrp-2);
322  double sn = sdat->value(nsc-1,nrp-1);
323  DbgLv(1) << "w:" << my_rank << ": s(01hmn)" << s0 << s1 << sh << sm << sn;
324  if ( dataset_count > 1 ) {
325  edat = &data_sets[offset+1]->run_data;
326  int nxx = nsc;
327  nsc = edat->scanCount();
328  nrp = edat->pointCount();
329  d0 = edat->scanData[0].rvalues[0];
330  d1 = edat->scanData[0].rvalues[1];
331  dh = edat->scanData[nsc/2].rvalues[nrp/2];
332  dm = edat->scanData[nsc-1].rvalues[nrp-2];
333  dn = edat->scanData[nsc-1].rvalues[nrp-1];
334  DbgLv(1) << "w:" << my_rank << ":d2(01hmn)" << d0 << d1 << dh << dm << dn;
335  dt = 0.0;
336  for ( int ss=0;ss<nsc;ss++ )
337  for ( int rr=0;rr<nrp;rr++ ) dt += edat->scanData[ss].rvalues[rr];
338  DbgLv(1) << "w:" << my_rank << ":dtot" << dt;
339  s0 = sdat->value(nxx+0,0);
340  s1 = sdat->value(nxx+0,1);
341  sh = sdat->value(nxx+nsc/2,nrp/2);
342  sm = sdat->value(nxx+nsc-1,nrp-2);
343  sn = sdat->value(nxx+nsc-1,nrp-1);
344  DbgLv(1) << "w:" << my_rank << ": s2(01hmn)" << s0 << s1 << sh << sm << sn;
345  }
346  DbgLv(1) << "w:" << my_rank << ": nsoli nsolo" << nsoli << nsolo;
347  DbgLv(1) << "w:" << my_rank << ": simpt men bott temp coef1"
348  << dset->simparams.simpoints
349  << dset->simparams.meniscus
350  << dset->simparams.bottom
351  << dset->simparams.temperature
352  << dset->simparams.rotorcoeffs[ 0 ];
353  DbgLv(1) << "w:" << my_rank << ": vbar soltype manual visc dens"
354  << dset->vbar20 << dset->solute_type << dset->manual
355  << dset->viscosity << dset->density;
356  DbgLv(1) << "w:" << my_rank << ": noisf alpha" << simu_values.noisflag
357  << simu_values.alpha << "s20w_c D20w_c vbar" << dset->s20w_correction
358  << dset->D20w_correction << dset->vbar20;
359  if ( dataset_count > 1 ) {
360  dset = data_sets[offset+1];
361  DbgLv(1) << "w:" << my_rank << ": 2)simpt men bott temp coef1"
362  << dset->simparams.simpoints
363  << dset->simparams.meniscus
364  << dset->simparams.bottom
365  << dset->simparams.temperature
366  << dset->simparams.rotorcoeffs[ 0 ];
367  DbgLv(1) << "w:" << my_rank << ": 2)vbar soltype manual visc dens"
368  << dset->vbar20 << dset->solute_type << dset->manual
369  << dset->viscosity << dset->density;
370  DbgLv(1) << "w:" << my_rank << ": 2)noisf alpha" << simu_values.noisflag
371  << simu_values.alpha << "s20w_c D20w_c vbar" << dset->s20w_correction
372  << dset->D20w_correction << dset->vbar20;
373  }
374  int nn = isols.size() - 1;
375  int mm = nn/2;
376  DbgLv(1) << "w:" << my_rank << ": sol0 solm soln" << isols[0].s << isols[0].k
377  << isols[mm].s << isols[mm].k << isols[nn].s << isols[nn].k;
378 }
379 //*DEBUG*
380 
381 }
382 
Generated on Sun Jun 21 2015 11:40:40 for UltraScan III by   doxygen 1.8.3.1-20130324