source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareParmela.cc @ 472

Last change on this file since 472 was 472, checked in by garnier, 10 years ago

Modification pour remettre en marche le Run. Desormais Transport passe, mais aucun test ne permet de dire si c est bon

File size: 19.6 KB
Line 
1#include "softwareParmela.h"
2#include "abstractElement.h"
3#include "parmelaParticle.h"
4#include "mathematicalConstants.h"
5#include "PhysicalConstants.h"
6#include "dataManager.h"
7#include "mixedTools.h"
8
9softwareParmela::softwareParmela() : abstractSoftware()
10{
11  nameOfSoftware_ = nomDeLogiciel("parmela");
12}
13
14softwareParmela::softwareParmela(string inputFileName,sectionToExecute* sect, dataManager* data) : abstractSoftware(inputFileName, sect, data)
15{
16  nameOfSoftware_ = nomDeLogiciel("parmela");
17  registerElement(nomdElements::RFgun,TBoolOk);
18  registerElement(nomdElements::drift,TBoolOk);
19  registerElement(nomdElements::cell,TBoolOk);
20  registerElement(nomdElements::bend,TBoolOk);
21  registerElement(nomdElements::soleno,TBoolOk);
22  registerElement(nomdElements::fit,TBoolIgnore);
23  registerElement(nomdElements::snapshot,TBoolIgnore);
24}
25
26void softwareParmela::setRelativeParmelaElementIndices() {
27  relativeParmelaElementIndices_.clear();
28  relativeParmelaElementIndices_.resize(numeroFin_deprecated_ - numeroDeb_deprecated_ + 1, -1);
29  cout << " setRelativeParmelaElementIndices() taille a priori : " << relativeParmelaElementIndices_.size() << endl;
30  abstractElement* elPtr = getSectionToExecute()->getElements().front();
31
32  bool there_is_rfgun = ( elPtr->getNomdElement().getElementType() == nomdElements::RFgun );
33  unsigned offsetNumElem;
34  // les elements de parmela sont indexes de 1 à max, s'il n'y a pas de rfgun
35  if ( !there_is_rfgun ) {
36    offsetNumElem = numeroDeb_deprecated_ -1;
37    // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
38  } else {
39    offsetNumElem = numeroDeb_deprecated_;
40  }
41
42  // index du premier element de parmela
43  int id= numeroDeb_deprecated_ - offsetNumElem;
44  unsigned k;
45  unsigned curseur = 0;
46  for ( k=0; k < getSectionToExecute()->getElements().size() ; k++ ) {
47    abstractElement* elem = getSectionToExecute()->getElements()[k];
48    cout << " liste PARMELA no absolu " << k << " relatif provisoire " << relativeParmelaElementIndices_.at(curseur) << endl;
49    // if ( elem->is_accepted_by_software(nameOfSoftware_) == TBoolOk ) {
50    if ( doAcceptElement(elem->getNomdElement().getElementType() )  == TBoolOk ) {
51      relativeParmelaElementIndices_.at(curseur) = id;
52      cout << " mis a " << id << endl;
53      id++;
54    }
55    curseur++;
56  }
57}
58
59
60bool softwareParmela::createInputFile(particleBeam* beamBefore, string workingDir)
61{
62  unsigned int k;
63
64  setRelativeParmelaElementIndices();
65  string name = workingDir + inputFileName_;
66  ofstream outfile;
67  outfile.open(name.c_str(), ios::out);
68  if (!outfile) {
69    dataManager_->consoleMessage(" softwareParmela::createInputFile : error opening output stream " );
70    cerr << " softwareParmela::createInputFile : error opening output stream " << name << endl;
71    return false;
72  }
73 
74  abstractElement* elPtr = getSectionToExecute()->getElements().front();
75  double initalKineticEnergy = 0.0;
76  bool there_is_rfgun = (elPtr->getNomdElement().getElementType() == nomdElements::RFgun );
77 
78  if ( !there_is_rfgun ) {
79    string nameOfInput = workingDir + "parin.input0";
80    if ( !beamToParmela(nameOfInput,beamBefore) ) return false;
81    initalKineticEnergy = beamBefore->referenceKineticEnergyMeV();
82    // les elements de parmela sont indexes de 1 à max, s'il n'y a pas de rfgun
83    //    offsetNumElem_ = numeroDeb_ -1;
84  } else {
85    elPtr->setPhaseStep( dataManager_->getGlobalParameters()->getIntegrationStep() );
86    initalKineticEnergy = elPtr->getInitialKineticEnergy();
87    // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
88  }
89 
90  outfile << "TITLE" << endl;
91  outfile << " titre provisoire " << endl;
92  outfile << "RUN /n0=1 /ip=999 /freq=" << dataManager_->getGlobalParameters()->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
93  outfile << "OUTPUT 0" << endl;
94  if ( there_is_rfgun ) {
95    //    outfile << dataManager_->getElementPointerFromNumero(numeroDeb_deprecated_)->parmelaOutputFlow();
96    outfile << elementsData(elPtr->parametersToSoftware());
97  } else {
98    outfile << "INPUT 0 /NP=" << beamBefore->getNbParticles() << endl;
99  }
100 
101 // retrouver le sector !!
102  for ( k =1; k <= getSectionToExecute()->getElements().size(); k++)
103    {
104      outfile << elementsData(getSectionToExecute()->getElements()[k]->parametersToSoftware());
105    }
106   
107  outfile << "ZOUT" << endl;
108  outfile << "START /wt=0.0 /dwt=" << dataManager_->getGlobalParameters()->getIntegrationStep() << "  /nsteps=" << dataManager_->getGlobalParameters()->getNbSteps() << "  nsc=" << dataManager_->getGlobalParameters()->getScPeriod() << "  /nout=10" << endl;
109  outfile << "END" << endl;
110  outfile.close();
111  dataManager_->consoleMessage("fichier input termine pour PARMELA");
112  return true;
113}
114
115bool softwareParmela::execute(string workingDir)
116{
117  bool ExecuteStatus = true;
118 
119  ostringstream sortie;
120  sortie << " EXECUTION DE PARMELA " <<  endl;
121   
122  string parmelaJob = workingDir + "parmela";
123  parmelaJob += string("   ");
124  parmelaJob += workingDir;
125 
126  string resultOfRun;
127  bool success = launchJob(parmelaJob,resultOfRun);
128
129
130  //  sortie << resultOfRun << endl;
131  if ( !success ) {
132    sortie << " launching of parmela failed " << endl;
133    ExecuteStatus = false;
134  } else {
135    sortie << " successful launching of parmela " << endl;
136    cout << " execution parmela MARCHE " << endl;
137    string::size_type nn = (resultOfRun).find("normal");
138    if ( nn == string::npos )
139      {
140        sortie << " abnormal exit of parmela " << endl;
141        ExecuteStatus = false;
142      } 
143  }
144 
145  dataManager_->consoleMessage(sortie.str());
146  return ExecuteStatus;
147}
148
149bool softwareParmela::buildBeamAfterElements( string workingDir)
150{
151  bool result = true;
152
153  if ( !ComputationLimitsOk_deprecated() ) return false;
154  unsigned curseur;
155  for ( unsigned int k=0; k < getSectionToExecute()->getElements().size() ; k++ ) {
156    abstractElement* elem = getSectionToExecute()->getElements()[k];
157      if ( elem == NULL ) {
158        dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : null pointer on element " );
159        return  false;       
160      }
161
162      curseur = k - numeroDeb_deprecated_;
163
164      if ( relativeParmelaElementIndices_.at(curseur) < 0 ) {
165
166        // si l'element doit etre ignore de parmela, on renvoie sur le diag precedent
167        particleBeam* lastDiag = dataManager_->updateCurrentDiagnostic(false);
168
169        // if(elem->getNomdElement().getElementType() == snapshot) {
170        //   // si cet element est un snapshot, on organise la sortie correspondante
171        //   string* param = elem->getParametersString();
172        //   string cliche = workingDir + param[2].c_str();
173        //   if( beamToParmela(cliche,lastDiag) ) {
174        //       //         cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam["<<avantDernier<<"]"<<endl;
175        //     cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam[ ]"<<endl;
176        //   }
177        // }
178        // si le numero est reconnu de parmela
179      } else {
180         
181        // on initialise une nouvelle sortie diagnostic
182          particleBeam* newDiag = dataManager_->updateCurrentDiagnostic(true);
183          vector<double> centroid;
184          bareParticle refPart;
185          vector<bareParticle> particles;
186          double freq= dataManager_->getGlobalParameters()->getFrequency();
187          unsigned numeroParmel;
188          numeroParmel = (unsigned)relativeParmelaElementIndices_.at(curseur);
189          cout << " lecture PARMELA el no absolu " << k << " numero relatif " << numeroParmel << " nom " << elem->getNomdElement().getExpandedName() << endl;
190        // lecture sortie parmela
191        if(!beamFromParmela(workingDir,numeroParmel,freq,centroid,refPart,particles))
192          {
193            // si echec, fin
194            dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : failure in reading parmela result for element " + elem->getLabel() + " for unknown reason " );
195            return  false;
196          } else {
197          // si succes, on complete le diagnostic
198          newDiag->setWithParticles(centroid,refPart,particles);
199        }
200      }
201    }
202  return result;
203}
204bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroParmel, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles )
205{   
206  string nomfilefais = workingDir + "parmdesz";
207  cout << " nom fichier desz : " << nomfilefais << endl;
208  FILE *filefais = fopen(nomfilefais.c_str(), "r");
209 
210  if ( filefais == (FILE*)0 ) {
211    dataManager_->consoleMessage(" beamFromParmela() erreur a l'ouverture du fichier 'parmdesz'");
212    cerr << " beamFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;
213    return false;
214  } else
215    cout << " beamFromParmela() : ouverture du fichier " << nomfilefais << endl;
216 
217  parmelaParticle partic;
218  std::vector<parmelaParticle> faisceau;
219  int testNombrePartRef =0;
220  double phaseRef = 0.0;
221 
222  while( partic.readFromParmelaFile(filefais) > 0 )
223    {
224     
225      if ( partic.ne == (int)numeroParmel ) {
226        if ( partic.np == 1 ) {
227          // en principe on est sur la particule de reference
228          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
229            printf(" ATTENTION part. reference douteuse  \n");
230            partic.imprim();
231          }
232          phaseRef = partic.phi;
233         
234          // le 'z' est 'absolu' (le long de la trajectoire)
235          TRIDVECTOR  posRef(partic.xx,partic.yy,partic.z);
236          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
237          refPart = bareParticle(posRef, betagammaRef);
238          testNombrePartRef++;
239          if ( testNombrePartRef != 1 ) {
240            dataManager_->consoleMessage(" beamFromParmela : nombre de part. de ref different de 1 :");
241            cerr << " nombre de part. de ref different de 1 : " << testNombrePartRef << " !! " << endl;
242            return false;
243          }
244        }
245        faisceau.push_back(partic);
246      }
247    } //while
248 
249  if ( faisceau.size() == 0)
250    {
251      string stringNum = mixedTools::intToString( (int)numeroParmel );
252      dataManager_->consoleMessage("beamFromParmela echec lecture  element numero relatif  parmela : " + stringNum);
253      cerr << " beamFromParmela echec lecture  element numero relatif  parmela " << numeroParmel << endl;
254      return false;
255    }
256 
257  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
258  // avec freq en Mhz et dphi en degres et résultat en cm:
259  double FACTEUR =  83.3333;  // ameliorer la precision
260 
261  // pour l'instant on choisit un centroid nul;
262  centroid.clear();
263  centroid = vector<double>(6,0.0);
264 
265  particles.clear();
266  particles.resize(faisceau.size(), bareParticle());
267  //  double x,xp,y,yp;
268  double xp, yp;
269  double betagammaz;
270  //  double betaz;
271  double cdt;
272  double dephas;
273  double g;
274  TRIDVECTOR  pos;
275  TRIDVECTOR betagamma;
276  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
277  // sont donnes en radians
278  for (unsigned k = 0; k < faisceau.size(); k++) {
279    //    x= faisceau.at(k).xx;
280    xp=faisceau.at(k).xxp;
281    //    y= faisceau.at(k).yy;
282    yp=faisceau.at(k).yyp;
283    // dephasage par rapport a la reference
284    dephas = faisceau.at(k).phi - phaseRef; // degrés
285    g = faisceau.at(k).wz/EREST_MeV;
286    betagammaz = faisceau.at(k).begamz;
287    //    betaz = betagammaz/(g+1.0);
288    //    deltaz = FACTEUR * betaz * dephas / referencefrequency;
289    cdt = FACTEUR * dephas / referencefrequency;
290    //    x += xp * deltaz;
291    //    y += yp * deltaz;
292    pos.setComponents(faisceau.at(k).xx,faisceau.at(k).yy,cdt);
293    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
294    particles.at(k) = bareParticle(pos,betagamma);
295  }
296 
297  return true;
298}
299
300// sauvegarde d'un 'particleBeam' sur un fichier parmela, en guise d'INPUT
301// pour l'instant de nom standard 'parin.input0'
302bool softwareParmela::beamToParmela(string nameOfFile,particleBeam* beam)
303{
304    if ( !beam->particleRepresentationOk() ) {
305        dataManager_->consoleMessage("softwareParmela::beamToParmela : beam not in particles form : not yet programmed");
306        cout << " softwareParmela::beamToParmela : beam not in particles form : not yet programmed " << endl;
307        return false;
308    }
309   
310    ofstream outfile;
311    outfile.open(nameOfFile.c_str(),ios::out);
312    if (!outfile) {
313        dataManager_->consoleMessage(" softwareParmela::beamToParmela : error opening output stream ");
314        cerr << " softwareParmela::beamToParmela : error opening output stream " << nameOfFile << endl;
315        return false;
316    }
317   
318    //    const vector<bareParticle>& partic = beam->getParticleVector();
319    double weight = 1.0;
320    double xx,yy,zz;
321    double begamx, begamy, begamz;
322    //    TRIDVECTOR pos, begam;
323    double zmin = GRAND;
324    double zmax = -zmin; 
325    double cdtmin, cdtmax;
326    beam->ZrangeCdt(cdtmin, cdtmax);
327    cout << " softwareParmela::beamToParmela cdtmin = " << cdtmin << " cdtmax = " << cdtmax << endl;
328 
329    for (unsigned k = 0; k < beam->getNbParticles(); k++) {
330        // partic.at(k).getPosition().getComponents(xx,yy,zz);
331        // partic.at(k).getBetaGamma().getComponents(begamx,begamy,begamz);
332      beam->coordonneesDeployees(k, -cdtmax).getComponents(xx,yy,zz);
333      beam->betaGamma(k).getComponents(begamx,begamy,begamz);
334      if ( zz > zmax) zmax = zz;
335      if ( zz < zmin ) zmin = zz;
336        outfile << xx << " " << begamx << " " <<  yy << " " << begamy << " " << zz << " " << begamz  << " " << weight << endl;
337    }
338    outfile.close();
339    cout << " softwareParmela::beamToParmela zmn = " << zmin << " zmax = " << zmax << endl;
340    return true;
341}
342
343
344string softwareParmela::elementsData(const vector< pair<string, vector<string> > >& donnees) const {
345  unsigned k;
346  cout << " PASSAGE softwareParmela::elementsData " << endl;
347  if ( donnees.at(0).first != "labelsGenericSpecific" ) {
348    cout << " softwareParmela::elementsData ERROR : element badly defined " << endl;
349    return string();
350  }
351  string genericName = donnees.at(0).second.at(0);
352  if ( genericName == "rfgun" ) return rfgunData(donnees);
353  if ( genericName == "cell" ) return cellData(donnees);
354  if ( genericName == "drift" ) return driftData(donnees);
355  if ( genericName == "solnd" ) return solenoData(donnees);
356  if ( genericName == "bend" ) return bendData(donnees);
357  return string();
358}
359
360string softwareParmela::rfgunData(const vector< pair<string, vector<string> > >& donnees) const {
361
362  cout << " PASSAGE softwareParmela::rfgunData " << endl;
363  string nmacrop = "0";
364  string sigma_t = "0.0";
365  string sigma_r = "0.0";
366  string E_cin = "0.0";
367  string sigma_E = "0.0";
368  string phaseStep = "0.0";
369  string totalCharge = "0.0";
370
371  unsigned k;
372  for ( k=1; k < donnees.size(); k++) {
373    if ( donnees.at(k).first == "nbMacroparticles" ) { 
374      nmacrop = donnees.at(k).second.at(0);
375    } else if ( donnees.at(k).first == "sigmasTR" ) {
376      sigma_t = donnees.at(k).second.at(0);
377      sigma_r = donnees.at(k).second.at(1);
378    } else if ( donnees.at(k).first == "kineticE" ) {
379      E_cin = donnees.at(k).second.at(0);
380      sigma_E = donnees.at(k).second.at(1);
381    } else if ( donnees.at(k).first == "phaseStep" ) {
382      phaseStep = donnees.at(k).second.at(0);
383    } else if ( donnees.at(k).first == "totalCharge" ) {
384      totalCharge = donnees.at(k).second.at(0);
385    }
386  }
387  ostringstream sortie;
388  // on prend les troncatures tmax et rmax à 3 sigmas
389  sortie << "INPUT 10 /np=" << nmacrop << "  /sigt=" << sigma_t << endl;
390  sortie << "  /tmax=" << 3.3*atof(sigma_t.c_str()) << " /sigr=" << sigma_r << endl;
391  sortie << " /rmax=" << 3.0*atof(sigma_r.c_str()) << " /W0=" << E_cin << " /dw0=" << sigma_E << endl;
392  sortie << " /dwt=" << phaseStep << " /ran=2" << endl;
393
394  // on doit entrer le nb vrai de part. (avec signe moins)
395  sortie << "SCHEFF /beami=" << -fabs( atof(totalCharge.c_str()) )/ELECTRONANOCOULOMB << " /nprog=2 /point=-1.7"  << endl;
396 
397  return sortie.str();
398}
399
400string softwareParmela::cellData(const vector< pair<string, vector<string> > >& donnees) const {
401  cout << " PASSAGE softwareParmela::cellData " << endl;
402
403  string lenght = "0.0";
404  string aperture = "0.0";
405  string initialPhase = "0.0";
406  string acceleratingField = "0.0";
407  string phaseStepMax = "0.0";
408  string acceleratingShapeFile = "";
409  string focusingMagFile = "";
410  string offsetMag = "0.0";
411  string scaleFactor = "0.0";
412
413  unsigned k;
414  for ( k=1; k < donnees.size(); k++) {
415    if ( donnees.at(k).first == "lenghtAperture" ) { 
416      lenght = donnees.at(k).second.at(0);
417      aperture = donnees.at(k).second.at(1);
418    } else if ( donnees.at(k).first == "phaseInitialStepmax" ) {
419      initialPhase = donnees.at(k).second.at(0);
420      phaseStepMax = donnees.at(k).second.at(1);
421    } else if ( donnees.at(k).first == "fieldValueFile" ) {
422      acceleratingField = donnees.at(k).second.at(0);
423      acceleratingShapeFile = donnees.at(k).second.at(1);
424    } else if ( donnees.at(k).first == "MagFocusingFileOffsetScale") {
425      focusingMagFile = donnees.at(k).second.at(0);
426      offsetMag = donnees.at(k).second.at(1);
427      scaleFactor = donnees.at(k).second.at(2);
428    }
429  }
430  ostringstream sortie;
431    sortie << "CELL /l=" << lenght << "  /aper=" << aperture << endl;
432    sortie << "  /iout=1  /phi0=" << initialPhase << " /E0=" << acceleratingField << endl;
433    sortie << " /nc=1 /dwtmax=" << phaseStepMax << " /sym=-1" << endl;
434    sortie << "CFIELD 1" << endl;
435    sortie << acceleratingShapeFile << endl;
436    if ( focusingMagFile != "") {
437      sortie << "POISSON /zoff=" << offsetMag << " /rmult=" << scaleFactor << endl;
438      sortie << focusingMagFile << endl;
439    }
440  return sortie.str();
441
442}
443
444string softwareParmela::driftData(const vector< pair<string, vector<string> > >& donnees) const {
445  cout << " PASSAGE softwareParmela::driftData " << endl;
446
447  string lenght = "0.0";
448  string aperture = "0.0";
449
450  unsigned k;
451  for ( k=1; k < donnees.size(); k++) {
452    if ( donnees.at(k).first == "lenghtAperture" ) { 
453      lenght = donnees.at(k).second.at(0);
454      aperture = donnees.at(k).second.at(1);
455    }
456  }
457  ostringstream sortie;
458  sortie << "DRIFT /l=" << lenght << "  /aper=" << aperture << "  /iout=1" << endl;
459  return sortie.str();
460}
461
462string softwareParmela::solenoData(const vector< pair<string, vector<string> > >& donnees) const {
463  cout << " PASSAGE softwareParmela::solenoData " << endl;
464  string lenght = "0.0";
465  string aperture = "0.0";
466  string B0 = "0.0";
467
468  unsigned k;
469  for ( k=1; k < donnees.size(); k++) {
470    if ( donnees.at(k).first == "lenghtAperture" ) { 
471      lenght = donnees.at(k).second.at(0);
472      aperture = donnees.at(k).second.at(1);
473    } else if ( donnees.at(k).first == "field" ) {
474      B0 = donnees.at(k).second.at(0);
475    }
476  }
477    ostringstream sortie;
478    // on passe l'induction magnetique en Gauss
479    sortie << "SOLENOID /l=" << lenght << "  /aper=" << aperture << "  /iout=1 /h=" << 1000.*atof(B0.c_str())  << endl;
480    return sortie.str();
481
482}
483
484string softwareParmela::bendData(const vector< pair<string, vector<string> > >& donnees) const {
485  cout << " PASSAGE softwareParmela::bendData " << endl;
486
487  string lenght = "0.0";
488  string aperture = "0.0";
489  string angleDeg = "0.0";
490  string momentum = "0.0";
491  string beta1 = "0.0";
492  string beta2 = "0.0";
493  unsigned k;
494  for ( k=1; k < donnees.size(); k++) {
495    if ( donnees.at(k).first == "lenghtAperture" ) { 
496      lenght = donnees.at(k).second.at(0);
497      aperture = donnees.at(k).second.at(1);
498    } else if ( donnees.at(k).first == "angleDegre" ) {
499      angleDeg = donnees.at(k).second.at(0);
500    } else if ( donnees.at(k).first == "momentum" ) {
501      momentum = donnees.at(k).second.at(0);
502     } else if ( donnees.at(k).first == "rotatedFaces" ) {
503      beta1 = donnees.at(k).second.at(0);
504      beta2 = donnees.at(k).second.at(1);
505    }   
506  }
507
508  double ecin = atof(momentum.c_str())/EREST_MeV;
509    ecin = ecin*ecin + 1.;
510    ecin = EREST_MeV*(sqrt(ecin) - 1.0);
511
512    ostringstream sortie;
513    // il faut entrer l'energie cinetique
514    sortie << "BEND /l=" << lenght << "  / aper=" << aperture << "  / iout=1 / wr=" << ecin << " /alpha=" << angleDeg << " / beta1=" << beta1;
515    sortie << " / beta2=" << beta2  << endl;
516
517    return sortie.str();
518
519}
Note: See TracBrowser for help on using the repository browser.