00001
00005 #include <libxml/xmlmemory.h>
00006 #include <libxml/parser.h>
00007 #include <zlib.h>
00008 #include <math.h>
00009 #include <string.h>
00010 #include <iostream>
00011 using namespace std;
00012 #include <limits>
00013 #include <omp.h>
00014 #include "def.h"
00015 #include "io.h"
00016 #include "run.h"
00017 #include "list.h"
00018 #include "collection.h"
00019 #include "grid.h"
00020 #include "container.h"
00021 #include "model.h"
00022 #include "species.h"
00023 using namespace Species;
00024 #include "domain.h"
00025 #include "gui.h"
00026
00027 void Boundary::init(int i, int j, REAL ymin[], REAL ymax[], int nspecies) {
00028 idir=i; iside=j; xmin=ymin; xmax=ymax;
00029 reactions=new Reaction[nspecies];
00030 }
00032 void Boundary::Inject(
00033 int ispecie,
00034 REAL vavx,
00035 REAL temperature,
00036 Collection<Molecule> *molecules
00037 ) {
00038 Molecule *molecule=molecules->append();
00039 molecule->Type(ispecie);
00040
00041
00042 REAL velx=vavx*fabs(GAUSS);
00043
00044 REAL drift=numeric_limits<REAL>::epsilon()+velx*Run::time.step;
00045 if(iside==0) {
00046 molecule->Coordinate(idir,xmin[idir]+drift);
00047 molecule->Velocity(idir,velx);
00048 } else if(iside==1) {
00049 molecule->Coordinate(idir,xmax[idir]-drift);
00050 molecule->Velocity(idir,-velx);
00051 } else {
00052 cerr<<"Inject: boundary not initialized\n";
00053 exit(1);
00054 }
00055 molecule->Time(Run::time.current);
00056 REAL dx=0.5*species[ispecie].Size();
00057 molecule->TimeStep(velx>SMALL?dx/vavx:0.0);
00058 for(int i=1;i<DIM;i++) {
00059 int j=(idir+i)%DIM;
00060 REAL
00061 x=xmin[j]+RND*(xmax[j]-xmin[j]),
00062
00063
00064 v=vavx*GAUSS;
00065 molecule->Coordinate(j,x);
00066 molecule->Velocity(j,v);
00067 }
00068 molecule->Temperature(temperature);
00069 GridContainer::put(molecule);
00070 }
00072 Domain::Domain(char *filename) {
00074 using namespace IO;
00075 using namespace Species;
00076 char *doctype=(char*)DOCTYPE;
00077 char *inpfilename=Run::configfile;
00079 char *p,word[MAXLINLEN];
00080 FILE *inp;
00081 for(int i=0;i<maxBoundaryTypes;i++)
00082 strcpy(boundaryName[i],"N/A");
00083 if(Run::option.verbose)printf("Domain::inpfilename=%s\n",inpfilename);
00084 xmlDocPtr doc;
00085 xmlNodePtr root;
00086 if (Run::option.verbose|Run::option.debug)
00087 { printf
00088 ( "Reading %s\n",
00089 inpfilename
00090 );fflush(stdout);
00091 }
00092 doc = xmlParseFile(inpfilename);
00093 if (doc == NULL )
00094 { fprintf(stderr,"XML parser failed in %s\n",inpfilename);
00095 exit(1);
00096 }
00097 root = xmlDocGetRootElement(doc);
00098 if (root == NULL)
00099 { fprintf(stderr,"Empty document: %s\n",inpfilename);
00100 xmlFreeDoc(doc);
00101 return;
00102 }
00103 if (xmlStrcmp(root->name, (const xmlChar *) doctype))
00104 { fprintf
00105 ( stderr,
00106 "document of the wrong type, root node != %s in %s\n",
00107 doctype,inpfilename
00108 );
00109 xmlFreeDoc(doc);
00110 return;
00111 }
00113 for
00114 ( xmlNodePtr cur = root->xmlChildrenNode;
00115 cur != NULL; cur = cur->next
00116 ) {
00118 if ((!xmlStrcmp(cur->name, (const xmlChar *)"xterm"))) {
00119 char str[WORDLENGTH];
00120 xmlChar *key=xmlNodeListGetString(doc,cur->xmlChildrenNode, 1);
00121 strncpy(str,(char*)key,WORDLENGTH);
00122 xmlFree(key);
00123 int xterm=0;
00124 sscanf(str,"%d",&xterm);
00125 if(xterm!=0)Run::option.xterm=true;
00126 else Run::option.xterm=false;
00127 if(Run::option.verbose) {
00128 printf("Xterm output: %d\n",xterm);fflush(stdout);
00129 }
00130 }
00132 if ((!xmlStrcmp(cur->name, (const xmlChar *)"molecules"))) {
00133 char str[WORDLENGTH];
00134 xmlChar *key=xmlNodeListGetString(doc,cur->xmlChildrenNode, 1);
00135 strncpy(str,(char*)key,WORDLENGTH);
00136 xmlFree(key);
00137 int nmolecules=0;
00138 sscanf(str,"%d",&nmolecules);
00139 if(Run::option.verbose) {
00140 printf("Maximum number of molecules: %d\n",nmolecules);fflush(stdout);
00141 }
00142 if(nmolecules>0) pool=new Pool<Molecule>(nmolecules);
00143 else pool=new Pool<Molecule>(1);
00144 molecules.init(nmolecules);
00145 }
00147 if ((!xmlStrcmp(cur->name, (const xmlChar *)"species"))) {
00148 char title[MAXLINLEN];
00149 if(getCharAttr(cur,(char*)"title",title)==0) {
00150 fprintf
00151 ( stderr,"CAN'T GET %s FOR %s IN %s\nAborting\n",
00152 "name","domain",inpfilename
00153 ); exit(1);
00154 }
00155 if(Run::option.verbose||Run::option.debug) printf("species: %s\n",title);
00156 xmlNodePtr tag=cur->xmlChildrenNode;
00157
00158 nspecies=0;
00159 while (tag != NULL) {
00160 if ((!xmlStrcmp(tag->name, (const xmlChar *)"specie")))
00161 nspecies++;
00162 tag=tag->next;
00163 }
00164 if(Run::option.verbose) printf("Number of species: %d\n",nspecies);
00165
00166 species=new Specie[nspecies+1];
00167 distribution=new int[nspecies];
00168
00169 tag=cur->xmlChildrenNode;
00170 for (int isp=0;tag != NULL&&isp<nspecies;tag = tag->next) {
00171 Specie *specie=species+isp;
00172 if ((!xmlStrcmp(tag->name, (const xmlChar *)"specie"))) {
00173 xmlChar *key=xmlGetProp(tag,(const xmlChar *)"id");
00174 if(key!=NULL) {
00175 char id[WORDLENGTH];
00176 strncpy(id,(char*)key,WORDLENGTH);
00177 strncpy(specie->Id(),id,WORDLENGTH);
00178 }
00179 REAL mass=parseFloat(doc,tag,(char*)"mass");
00180 specie->Mass((REAL)mass);
00181 REAL size=parseFloat(doc,tag,(char*)"size");
00182 specie->Size((REAL)size);
00183
00184
00185
00186
00187
00188
00189
00190
00191 REAL cp=parseFloat(doc,tag,(char*)"cp");
00192 if(cp<=2.5*GasConstant)cp=2.5*GasConstant;
00193 specie->Cp(cp);
00194 if(Run::option.verbose){
00195 printf(
00196 "specie %d: name=%s, mass='%g', size='%g', cp=%g\n",
00197 isp,species[isp].Id(),species[isp].Mass(),species[isp].Size(),cp
00198 );fflush(stdout);
00199 }
00200 isp++;
00201 }
00202 }
00203 strcpy(species[VOIDSPECIE].Id(),"VOID\0");
00204 species[VOIDSPECIE].Size(TINY);
00205 species[VOIDSPECIE].Mass(0.0);
00206
00207 reactions=new Reaction[nspecies*nspecies];
00209 while(tag!=NULL) {
00210 if (!xmlStrcmp(tag->name, (const xmlChar *)"reaction")) {
00211 char reactants[WORDLENGTH],products[WORDLENGTH];
00212 xmlChar *key=xmlGetProp(tag,(const xmlChar *)"reactants");
00213 if(key==NULL){cerr<<"Keyword 'reactants' missing in reaction tag\n";exit(1);}
00214 strncpy(reactants,(char*)key,WORDLENGTH);
00215 key=xmlGetProp(tag,(const xmlChar *)"products");
00216 if(key==NULL){cerr<<"Keyword 'products' missing in reaction tag\n";exit(1);}
00217 strncpy(products,(char*)key,WORDLENGTH);
00218 REAL
00219 activationEnergy=BoltzmannConstant/AtomicMassUnit*parseFloat(doc,tag,(char*)"activation"),
00220 probability=parseFloat(doc,tag,(char*)"probability"),
00221
00222 enthalpy=1.0e3*parseFloat(doc,tag,(char*)"enthalpy")/(AvogadroNumber*AtomicMassUnit);
00223
00224 char *r=reactants,*p=products,word[WORDLENGTH],*w=word;
00225 while(isspace(*r))r++; while(isspace(*p))p++;
00226 int ir0=0,ir1=0,ip0=0,ip1=0,j=0;
00227 for(w=word;w-word<WORDLENGTH-1&&!isspace(*r);r++,w++) *w=*r; *w='\0';
00228 for(ir0=0;ir0<nspecies;ir0++) if(strncmp(species[ir0].Id(),word,WORDLENGTH)==0) break;
00229 if(ir0==nspecies){cerr<<"First reactant not found for '"<<word<<"'\n";exit(1);}
00230 while(isspace(*r))r++;for(w=word;w-word<WORDLENGTH-1&&!isspace(*r);r++,w++) *w=*r; *w='\0';
00231 for(ir1=0;ir1<nspecies;ir1++) if(strncmp(species[ir1].Id(),word,WORDLENGTH)==0) break;
00232 if(ir1==nspecies){cerr<<"Second reactant not found for '"<<word<<"'\n";exit(1);}
00233 while(isspace(*p))p++;for(w=word;w-word<WORDLENGTH-1&&!isspace(*p);p++,w++) *w=*p; *w='\0';
00234 for(ip0=0;ip0<nspecies;ip0++) if(strncmp(species[ip0].Id(),word,WORDLENGTH)==0) break;
00235 if(ip0==nspecies){cerr<<"Product not found for '"<<word<<"'\n";exit(1);}
00236 while(*p!='\0'&&isspace(*p))p++;
00237 if(*p!='\0') {
00238 for(w=word;w-word<WORDLENGTH-1&&!isspace(*p);p++,w++) *w=*p; *w='\0';
00239 for(ip1=0;ip1<nspecies;ip1++) if(strncmp(species[ip1].Id(),word,WORDLENGTH)==0) break;
00240 } else ip1=VOIDSPECIE;
00241 int ireact0=ir0*nspecies+ir1, ireact1=ir1*nspecies+ir0;
00242
00243 Reaction *react0=reactions+ireact0, *react1=reactions+ireact1;
00244 react0->Add(ip0,ip1,activationEnergy,probability,enthalpy);
00245 react1->Add(ip0,ip1,activationEnergy,probability,enthalpy);
00246
00247 }
00248 tag=tag->next;
00249 }
00250 if(Run::option.verbose||Run::option.debug) {
00251 cout<<"Reactions between "<<nspecies<<" species:\n";
00252 for(int i=0;i<nspecies;i++)
00253 for(int j=0;j<nspecies;j++){
00254 Reaction *r=reactions+i*nspecies+j;
00255 Reaction::Outcome *o=r->First();
00256 if(o!=NULL)
00257 do{
00258 int ip0=o->Product(0), ip1=o->Product(1);
00259 REAL rate=o->Probability(), enthalpy=o->Enthalpy();
00260 if(rate>0.0) {
00261 cout<<i<<','<<j<<':'<<species[i].Id()<<" + "<<species[j].Id()<<" = "
00262 <<species[ip0].Id()<<" + "<<species[ip1].Id()
00263 <<"; rate="<<rate<<", ent="<<enthalpy<<endl;
00264 cout.flush();
00265 }
00266 } while((o=r->Next())!=NULL);
00267 }
00268 }
00269 }
00271 if ((!xmlStrcmp(cur->name, (const xmlChar *)"domain"))) {
00272 char name[MAXLINLEN];
00273 if(getCharAttr(cur,(char*)"name",name)==0) {
00274 fprintf
00275 ( stderr,"CAN'T GET name FOR domain IN %s\nAborting\n",
00276 inpfilename
00277 ); exit(1);
00278 }
00279 if(Run::option.verbose) printf("Domain: %s\n",name);
00280 if(parseWord(doc,cur,(char*)"type",word)==0) {
00281 fprintf
00282 ( stderr,"Tag \"type\" is missing in %s\nAborting\n",
00283 inpfilename
00284 ); exit(1);
00285 }
00286 if(Run::option.verbose){printf("Domain type: %s\n",word);fflush(stdout);}
00287 if(strcmp(word,"box")!=0) {
00288 fprintf(stderr,"ERROR: Only domain type 'box' is supported\n");
00289 exit(1);
00290 }
00291
00292
00293 strcpy(boundaryName[insideBoundary],"inside");
00294 strcpy(boundaryName[periodicBoundary],"periodic");
00295 strcpy(boundaryName[elasticBoundary],"elastic");
00296 strcpy(boundaryName[openBoundary],"open");
00297 for
00298 ( xmlNodePtr next=cur->xmlChildrenNode;
00299 next != NULL;
00300 next=next->next
00301 ) {
00302 if ((!xmlStrcmp(next->name, (const xmlChar *)"grid"))) {
00303
00304 if(Run::option.verbose){printf("Particle partition grid:\n");fflush(stdout);}
00305 if(parseWord(doc,next,(char*)"cellsize",word)!=0) {
00306 REAL cellsize=0.0;
00307 sscanf(word,"%lg", &cellsize);
00308 GridContainer::cellsize=cellsize;
00309 if(Run::option.verbose) printf ( "\tiCell size: %g\n",cellsize);
00310 }
00311 } else
00312 if ((!xmlStrcmp(next->name, (const xmlChar *)"energy"))) {
00313 if(Run::option.verbose){printf("Interaction Potential:\n");fflush(stdout);}
00314 if(parseWord(doc,next,(char*)"lengthscale",word)!=0) {
00315 REAL lengthscale=0.0,strength=0.0;
00316 sscanf(word,"%lg", &lengthscale);
00317 Potential::lengthscale(lengthscale);
00318 if(Run::option.verbose) printf ( "\tlengthscale: %g\n",lengthscale);
00319 }
00320 if(parseWord(doc,next,(char*)"strength",word)!=0) {
00321 REAL strength=0.0;
00322 sscanf(word,"%lg", &strength);
00323 Potential::strength(strength);
00324 if(Run::option.verbose) printf("\tstrength: %g\n",strength);
00325 }
00326 }
00327 else
00328 if((!xmlStrcmp(next->name, (const xmlChar*)"bounds"))) {
00329 REAL *xmin=bounds[0],*xmax=bounds[1];
00330 xmlChar *key=xmlNodeListGetString(doc,next->xmlChildrenNode,1);
00331 strncpy(word,(char*)key,MAXLINLEN);
00332 xmlFree(key);
00333 sscanf
00334 ( word,"%lg %lg %lg %lg %lg %lg",
00335 xmin,xmax,xmin+1,xmax+1,xmin+2,xmax+2
00336 );
00337 if(Run::option.verbose)
00338 printf
00339 ( "\tbounds: (%lg:%lg), (%lg:%lg), (%lg:%lg)\n",
00340 xmin[0],xmax[0],xmin[1],xmax[1],xmin[2],xmax[2]
00341 );
00342
00343 bool adjusted=false;
00344 for(int i=0;i<DIM;i++) {
00345 if(xmin[i]>minBound(i)) {
00346 xmin[i]=minBound(i);
00347 adjusted=true;
00348 } else setMinBound(i,xmin[i]);
00349 if(xmax[i]<maxBound(i)) {
00350 xmax[i]=maxBound(i);
00351 adjusted=true;
00352 } else setMaxBound(i,xmax[i]);
00353 }
00354 if(adjusted)
00355 printf
00356 ( "\tDomain bounds adjusted: (%lg:%lg), (%lg:%lg), (%lg:%lg)\n",
00357 xmin[0],xmax[0],xmin[1],xmax[1],xmin[2],xmax[2]
00358 );
00359
00360 for(int idir=0;idir<DIM;idir++) {
00361 int i1=(idir+1)%DIM, i2=(idir+2)%DIM;
00362 REAL area=(xmax[i1]-xmin[i1])*(xmax[i2]-xmin[i2]);
00363 for(int iside=0;iside<2;iside++) {
00364 boundaries[idir][iside].area=area;
00365 if(Run::option.verbose)
00366 cout<<"Boundary ("<<idir<<','<<iside<<").area="<<area<<endl;
00367 }
00368 }
00369 }
00370 else
00371 if((!xmlStrcmp(next->name, (const xmlChar*)"bulk"))) {
00372 if(Run::option.verbose)printf("Domain bulk properties:\n");
00373 REAL temperature=300.0;
00374 if(parseFloat(doc,next,(char*)"temperature",temperature)!=0) {
00375 bulk.Temperature(temperature);
00376 if(Run::option.verbose)printf("\tbulk temperature=%g\n",temperature);
00377 }
00378 if(temperature<=0.0) {
00379 printf("Wrong domain temperature: %g\n",temperature);
00380 exit(1);
00381 }
00382
00383 xmlNodePtr tag=next->xmlChildrenNode;
00384 while(tag!=NULL) {
00385 if(!xmlStrcmp(tag->name,(const xmlChar*)"specie")) {
00386
00387 char id[MAXLINLEN];
00388 if(getCharAttr(tag,(char*)"id",id)) {
00389 if(Run::option.verbose||Run::option.debug) {
00390 cout<<"\tBulk Specie id="<<id<<endl;cout.flush();
00391 }
00392
00393 Specie *specie=NULL;
00394 for(int isp=0;isp<nspecies;isp++) {
00395 Specie *spc=species+isp;
00396 if(strcmp(id,spc->Id())==0){
00397 specie=spc;break;
00398 }
00399 }
00400 if(specie==NULL) {
00401 cerr<<"FAILED TO LOCATE SPECIE ID="<<id<<endl;
00402 exit(1);
00403 }
00404
00405 Gas *gas=new Gas(specie);
00406 bulk.gases.append(gas);
00407 REAL density=0.0;
00408 if(parseFloat(doc,tag,(char*)"density",density))
00409 if(Run::option.verbose||Run::option.debug) {
00410 cout<<"\tDensity [kg/m^3]="<<density<<endl;
00411 }
00412 else
00413 if(Run::option.verbose||Run::option.debug) {
00414 cout<<"Density not specified for specie "<<id<<", assumed 0\n";
00415 }
00416 gas->density=density;
00417 } else {
00418 cerr<<"No Specie ID specified in "<<inpfilename<<endl;
00419 exit(1);
00420 }
00421 }
00422 tag=tag->next;
00423 }
00424 }
00425 else
00426 if((!xmlStrcmp(next->name, (const xmlChar*)"boundary"))) {
00427 char id[MAXLINLEN];
00428 REAL *xmin=bounds[0],*xmax=bounds[1];
00429 if(getCharAttr(next,(char*)"id",id)==0) {
00430 fprintf
00431 ( stderr,"CAN'T GET id FOR boundary IN %s\nAborting\n",
00432 inpfilename
00433 ); exit(1);
00434 }
00435 int idir=-1,iside=-1;
00436 if(strcmp("left",id)==0) {idir=0;iside=0;}
00437 else if(strcmp("right",id)==0) {idir=0;iside=1;}
00438 else if(strcmp("bottom",id)==0) {idir=1;iside=0;}
00439 else if(strcmp("top",id)==0) {idir=1;iside=1;}
00440 else if(strcmp("front",id)==0) {idir=2;iside=0;}
00441 else if(strcmp("rear",id)==0) {idir=2;iside=1;}
00442 else {
00443 fprintf(stderr,"UNKNOWN BOUNDARY ID '%' in %s\n",id,inpfilename);
00444 exit(1);
00445 }
00446 Boundary *boundary=&boundaries[idir][iside];
00447 boundary->init(idir,iside,xmin,xmax,nspecies);
00448 if(parseWord(doc,next,(char*)"type",word)!=0) {
00449 int itype=0;
00450 p=word;
00451 while(isspace(*p))p++;
00452 for(itype=0;itype<maxBoundaryTypes;itype++) {
00453 char *name=boundaryName[itype];
00454 if(strncmp(p,name,strlen(name))==0) break;
00455 }
00456 if(itype==maxBoundaryTypes) {
00457 fprintf(stderr,"Can't assign boundary type '%s'\n",p);
00458 exit(1);
00459 }
00460 boundary->type=(BoundaryTypes)itype;
00461 }
00462 REAL temperature=-10.0;
00463 if(parseFloat(doc,next,(char*)"temperature",temperature)!=0) {
00464 boundary->temperature=temperature;
00465 }
00466 if(temperature<=0.0) boundary->adiabatic=true;
00467 else boundary->adiabatic=false;
00468 if(Run::option.verbose) {
00469 char *type=(char*)"isothermal";
00470 if(boundary->adiabatic) type=(char*)"adiabatic";
00471 printf(
00472 "Boundary %s (%d,%d): type=%d: %s, temperature=%g, heat-transfer: %s\n",
00473 id,idir,iside,(int)boundary->type,boundaryName[boundary->type],boundary->temperature,type
00474 );
00475 }
00476 xmlNodePtr tag=next->xmlChildrenNode;
00477 while(tag!=NULL) {
00478 if(!xmlStrcmp(tag->name,(const xmlChar*)"specie")) {
00479
00480 char id[MAXLINLEN];
00481 if(getCharAttr(tag,(char*)"id",id)) {
00482 if(Run::option.verbose||Run::option.debug) {
00483 cout<<"\tBoundary Specie id="<<id<<endl;cout.flush();
00484 }
00485
00486 Specie *specie=NULL;
00487 for(int isp=0;isp<nspecies;isp++) {
00488 Specie *spc=species+isp;
00489 if(strcmp(id,spc->Id())==0){
00490 specie=spc;break;
00491 }
00492 }
00493 if(specie==NULL) {
00494 cerr<<"FAILED TO LOCATE SPECIE ID="<<id<<endl;
00495 exit(1);
00496 }
00497
00498
00499 Gas *gas=new Gas(specie);
00500 boundary->gases.append(gas);
00501 REAL density=0.0;
00502 if(parseFloat(doc,tag,(char*)"density",density))
00503 if(Run::option.verbose||Run::option.debug) {
00504 cout<<"\tDensity [kg/m^3]="<<density<<endl;
00505 }
00506 else
00507 if(Run::option.verbose||Run::option.debug) {
00508 cout<<"Density not specified for specie "<<id<<", assumed 0\n";
00509 }
00510 gas->density=density;
00511 } else {
00512 cerr<<"No Specie ID specified in "<<inpfilename<<endl;
00513 exit(1);
00514 }
00515 }
00516 else if(!xmlStrcmp(tag->name,(const xmlChar*)"reaction")) {
00518
00519 char reactant[WORDLENGTH],products[WORDLENGTH];
00520 xmlChar *key=xmlGetProp(tag,(const xmlChar *)"reactant");
00521 if(key==NULL){cerr<<"Keyword 'reactant' missing in boundary reaction tag\n";exit(1);}
00522 strncpy(reactant,(char*)key,WORDLENGTH);
00523 if(Run::option.verbose||Run::option.debug) {
00524 cout<<"\tBoundary Reactant ="<<reactant<<endl;cout.flush();
00525 }
00526 key=xmlGetProp(tag,(const xmlChar *)"products");
00527 if(key==NULL){cerr<<"Keyword 'products' missing in boundary reaction tag\n";exit(1);}
00528 strncpy(products,(char*)key,WORDLENGTH);
00529 if(Run::option.verbose||Run::option.debug) {
00530 cout<<"\tBoundary Products ="<<products<<endl;cout.flush();
00531 }
00532 REAL
00533 activationEnergy=BoltzmannConstant/AtomicMassUnit*parseFloat(doc,tag,(char*)"activation"),
00534 probability=parseFloat(doc,tag,(char*)"probability"),
00538 #ifdef ADSORPTION
00539 residencetime=parseFloat(doc,tag,(char*)"time"),
00540 #endif
00541 enthalpy=parseFloat(doc,tag,(char*)"enthalpy");
00542
00543 char *r=reactant,*p=products,word[WORDLENGTH],*w=word;
00544 while(isspace(*r))r++; while(isspace(*p))p++;
00545 int ir0=0,ip0=0,ip1=0,j=0;
00546
00547 for(w=word;w-word<WORDLENGTH-1&&!isspace(*r);r++,w++) *w=*r; *w='\0';
00548 for(ir0=0;ir0<nspecies;ir0++) if(strncmp(species[ir0].Id(),word,WORDLENGTH)==0) break;
00549 if(ir0==nspecies){cerr<<"Boundary reactant not found for '"<<word<<"'\n";exit(1);}
00550
00551 while(isspace(*p))p++;for(w=word;w-word<WORDLENGTH-1&&!isspace(*p);p++,w++) *w=*p; *w='\0';
00552 for(ip0=0;ip0<nspecies;ip0++) if(strncmp(species[ip0].Id(),word,WORDLENGTH)==0) break;
00553 if(ip0==nspecies){cerr<<"Boundary products not found for '"<<word<<"'\n";exit(1);}
00554 while(*p!='\0'&&isspace(*p))p++;
00555 if(*p!='\0') {
00556 for(w=word;w-word<WORDLENGTH-1&&!isspace(*p);p++,w++) *w=*p; *w='\0';
00557 for(ip1=0;ip1<nspecies;ip1++) if(strncmp(species[ip1].Id(),word,WORDLENGTH)==0) break;
00558 } else ip1=VOIDSPECIE;
00559 #ifdef ADSORPTION
00560 boundary->reactions[ir0].Add(ip0,ip1,activationEnergy,probability,enthalpy,residencetime);
00561 #else
00562 boundary->reactions[ir0].Add(ip0,ip1,activationEnergy,probability,enthalpy);
00563 #endif
00564 }
00565 tag=tag->next;
00566 }
00567 if(Run::option.verbose||Run::option.debug) {
00568 cout<<"Reactions between "<<nspecies<<" boundary species:\n";
00569 for(int i=0;i<nspecies;i++){
00570 Reaction *r=boundary->reactions+i;
00571 Reaction::Outcome *o=r->First();
00572 if(o!=NULL) {
00573 do{
00574 int ip0=o->Product(0), ip1=o->Product(1);
00575 REAL rate=o->Probability(), enthalpy=o->Enthalpy();
00576 if(rate>0.0) {
00577 cout<<i<<':'<<species[i].Id()<<" = "
00578 <<species[ip0].Id()<<" + "<<species[ip1].Id()
00579 <<"; probability="<<rate<<", enthalpy="<<enthalpy<<endl;
00580 cout.flush();
00581 }
00582 } while((o=r->Next())!=NULL);
00583 }
00584 }
00585 }
00586 }
00587 }
00588 }
00589 }
00590 xmlFreeDoc(doc);
00591 if(Run::option.verbose){printf("File %s parsed\n",inpfilename);fflush(stdout);}
00592 bulk.init(XMIN,XMAX);
00593 if(Run::option.restart)
00594 load(filename);
00595 else
00596 init();
00597 #ifdef LOCAL
00598 GridContainer::init(
00599 XMIN, XMAX,
00600
00601 GridContainer::cellsize,
00602 pool
00603 );
00604 if(Run::option.verbose||Run::option.debug) {
00605 if(Run::option.debug)cout<<"pool:"<<pool<<endl;cout.flush();
00606 cout<<"GridContainer::init:molecules:"<<molecules.number()<<endl;cout.flush();
00607 }
00608 if(molecules.number()>0) {
00609 molecules.goFirst();
00610 do {
00611 GridContainer::put(molecules.Current());
00612 molecules.goNext();
00613 } while(!molecules.isFirst());
00614 }
00615 #endif
00616 }
00617 Domain::~Domain() {
00618 if(Run::option.debug) {cout<<"Deleting domain\n";cout.flush();}
00619
00620 }
00621 void Domain::init() {
00622 double temperature=bulk.Temperature(),
00623 volume=(double)bulk.Volume()*1e-27;
00624 List<Gas> *gases=&bulk.gases;
00625 int imolecule,nmolecules=0;
00626 if(gases->number()>0) {
00627 gases->goFirst();
00628 do {
00629 Gas *gas=gases->Current();
00630 Specie *specie=gas->specie;
00631
00634 double
00635 mu=(double)specie->Mass(),
00636 mass=(double)AtomicMassUnit*mu,
00637
00638
00639 energy=BoltzmannConstant*temperature,
00640 Vavx=sqrt(energy/mass),
00641 density=(double)gas->density,
00642 totmass=density*volume;
00643 int nmol=(int)round(totmass/mass);
00644
00645 if(Run::option.verbose) {
00646 printf("Domain::init:Injecting %d molecules of specie %s into the bulk\n",nmol,specie->Id());
00647 }
00648 for(int inj=0;inj<nmol;inj++) {
00649 Molecule *molecule = bulk.inject((int)(specie-Species::species),Vavx,&molecules);
00650 molecule->Temperature(temperature);
00651 }
00652 gases->goNext();
00653 } while(!gases->isFirst());
00654 }
00655 }
00657 void Bulk::init(REAL xmn[], REAL xmx[]) {
00658 volume=1.0;
00659 for(int i=0;i<DIM;i++) {
00660 xmin[i]=xmn[i];xmax[i]=xmx[i];
00661 volume*=xmax[i]-xmin[i];
00662 }
00663 }
00665 Molecule *Bulk::inject(
00666 int ispecie,
00667 REAL velx,
00668 Collection<Molecule> *molecules
00669 ) {
00670 Molecule *molecule=molecules->append();
00671 molecule->Type(ispecie);
00672 molecule->Time(Run::time.current);
00673 REAL radius=0.5*species[ispecie].Size();
00674 molecule->TimeStep(velx>SMALL?radius/velx:0.0);
00675
00676 for(int i=0;i<DIM;i++) {
00677
00678 REAL x=xmin[i]+RND*(xmax[i]-xmin[i]);
00679 molecule->Coordinate(i,x);
00680
00681 REAL
00682 v=velx*GAUSS;
00683 molecule->Velocity(i,v);
00684 }
00685 return molecule;
00686 }
00687 void Domain::load(char *filename) {
00688 uLong nbuf=MAXLINLEN*sizeof(int);
00689 Byte *buf=(Byte*)calloc((uInt)nbuf,1);
00690 gzFile gzinp;
00691 if((gzinp = gzopen(filename, "rb"))==NULL)
00692 { cerr << "CAN'T OPEN "<<filename<<endl;cerr.flush();
00693 exit(1);
00694 }
00695 if(Run::option.verbose||Run::option.debug)
00696 cout<<"Reading from "<<filename<<endl;
00697 int nmolecules=0;
00698 while(true) {
00699 char *p;
00700 gzgets(gzinp, (char*)buf,nbuf);
00701 if((char)buf[0]!='#')break;
00702 char *s=(char*)buf;
00703 if((p=strstr(s,"time"))!=NULL){
00704 p+=5;
00705 while(isspace(*p))p++;
00706 REAL starttime=0.0;
00707 sscanf(p,"%lg",&starttime);
00708 Run::time.current=starttime;
00709 if(Run::time.start>starttime) Run::time.start=starttime;
00710 if(Run::option.verbose) {cout<<"Start time="<<starttime<<endl;cout.flush();}
00711 }
00712 if((p=strstr(s,"output"))!=NULL){
00713 p+=6;
00714 while(isspace(*p))p++;
00715 int ioutput=0;
00716 sscanf(p,"%d",&ioutput);
00717 IO::ioutput=ioutput;
00718 if(Run::option.verbose) {cout<<"Output dataset="<<IO::ioutput<<endl;cout.flush();}
00719 }
00720 if((p=strstr(s,"molecules"))!=NULL){
00721 p+=9;
00722 while(isspace(*p))p++;
00723 int ioutput=0;
00724 sscanf(p,"%d",&nmolecules);
00725 if(Run::option.verbose) {cout<<"Number of molecules="<<nmolecules<<endl;cout.flush();}
00726 }
00727 }
00728 int imolecule;
00729 REAL ymin[DIM],ymax[DIM];
00730 if(nmolecules>0) {
00731 for(int i=0;i<DIM;i++) {
00732 ymin[i]=LARGE;
00733 ymax[i]=-LARGE;
00734 }
00735 } else {
00736 for(int i=0;i<DIM;i++)ymin[i]=ymax[i]=0.0;
00737 }
00738 REAL maxsize=0.0;
00739
00740
00741 for(imolecule=0;
00742 imolecule<nmolecules&&!gzeof(gzinp);
00743 imolecule++,
00744 gzgets(gzinp, (char*)buf, nbuf)
00745 ) {
00746 char *s=(char*)buf;
00747 if(*s=='#') {imolecule--; continue;}
00748 Molecule *molecule=molecules.append();
00749 char id[WORDLENGTH];
00750 sscanf(s,"%s",&id);
00751 int type=-1;
00752 for(int i=0;i<nspecies;i++) {
00753 if(strncmp(id,species[i].Id(),WORDLENGTH)==0) {
00754 type=i;
00755 break;
00756 }
00757 }
00758 if(type<0) {
00759 cout<<"Unknown specie identifier "<<id<<endl;
00760 exit(1);
00761 }
00762 molecule->Type(type);
00763
00764 REAL size=species[type].Size();
00765 if(maxsize<size)maxsize=size;
00766
00767
00768
00769 while(!isspace(*++s));
00770 while(isspace(*++s));
00771 for (int i=0;i<DIM;i++) {
00772 REAL x;
00773 sscanf(s,"%lg",&x);
00774 molecule->Coordinate(i,x);
00775 if(ymin[i]>x)ymin[i]=x;
00776 if(ymax[i]<x)ymax[i]=x;
00777 while(!isspace(*++s));
00778 while(isspace(*++s));
00779 }
00780 for (int i=0;i<DIM;i++) {
00781 REAL v;
00782 sscanf(s,"%lg",&v);
00783 molecule->Velocity(i,v);
00784 while(!isspace(*++s));
00785 while(isspace(*++s));
00786 }
00787 { REAL e;
00788 sscanf(s,"%lg",&e);
00789 molecule->InternalEnergy(e);
00790 while(!isspace(*++s));
00791 while(isspace(*++s));
00792 }
00793 { REAL t,dt;
00794 sscanf(s,"%lg",&t);
00795 molecule->Time(t);
00796 while(!isspace(*++s));
00797 while(isspace(*++s));
00798 sscanf(s,"%lg",&dt);
00799 molecule->TimeStep(dt);
00800 while(!isspace(*++s));
00801 while(isspace(*++s));
00802 }
00803 #ifdef DEBUG
00804 if(Run::option.debug) {
00805 cout<<imolecule<<": "<<id<<", x:";cout.flush();
00806 for(int i=0;i<DIM;i++) {
00807 cout<<" "<<molecule->Coordinate(i);
00808 }
00809 cout<<imolecule<<", v:";cout.flush();
00810 for(int i=0;i<DIM;i++) {
00811 cout<<" "<<molecule->Velocity(i);
00812 }
00813 cout<<endl;cout.flush();
00814 }
00815 #endif
00816 }
00817 if(imolecule!=nmolecules) {
00818 cerr<<"Number of molecules specified ("<<nmolecules
00819 <<") does not agree with actually found ("<<imolecule<<")\n";
00820 cerr.flush();exit(1);
00821 }
00822 if(Potential::lengthscale()>0.0) Potential::Cutoff(Potential::lengthscale());
00823 else if(maxsize>0.0)Potential::Cutoff(2.0*maxsize);
00824 else { cerr<<"BAD LENGTHSCALE\n";cerr.flush();exit(1);}
00825 if(Run::option.debug) cout << molecules.number()
00826 << " molecules read. Interaction Cutoff = "<<Potential::Cutoff()<<endl;
00827 gzclose(gzinp);
00828 }
00829 int Domain::computeBounds(REAL x0[], REAL x1[]) {
00830
00831 for(int i=0;i<DIM;i++) {
00832 x0[i]=LARGE;
00833 x1[i]=-LARGE;
00834 }
00835 molecules.goFirst();
00836 int nmolecules=0;
00837 do {
00838 Molecule *molecule=molecules.Current();
00839 REAL *x=molecule->Coordinates();
00840 for(int i=0;i<DIM;i++) {
00841 if(x0[i]>x[i])x0[i]=x[i];
00842 if(x1[i]<x[i])x1[i]=x[i];
00843 }
00844 molecules.goNext(); nmolecules++;
00845 } while(!molecules.isFirst());
00846 return nmolecules;
00847 }
00849 int Domain::run(int niter) {
00850 int
00851 maxcounters=molecules.maxCounters();
00852 Run::nthreads=omp_get_max_threads();
00853 if(Run::nthreads>maxcounters) Run::nthreads=maxcounters;
00854 omp_set_num_threads(Run::nthreads);
00855 int *distribution=new int[nspecies];
00856 REAL currenttime=Run::time.current,
00857 starttime=Run::time.start,endtime=Run::time.end;
00858 if (Run::option.debug)
00859 cout << "Running "<< niter << " iterations ...";
00860 time_t start_wtime=time((time_t) NULL);
00861
00862 if(molecules.number()>0) {
00863 molecules.goFirst();
00864 int imolecule=0;
00865 do {
00866 Molecule *molecule=molecules.Current();
00867 if(molecule->Type()!=VOIDSPECIE)molecule->index(imolecule++);
00868 molecules.goNext();
00869 } while(!molecules.isFirst());
00870 }
00871 REAL maxjump=0.5*GridContainer::cellsize;
00872 int iter=0;
00874 for(;iter<niter&¤ttime<=endtime;iter++) {
00875 Run::time.step=LARGE;
00876 #pragma omp parallel
00877 {
00878 int ithread = omp_get_thread_num();
00879
00880 molecules.goFirst(ithread);
00881 do {
00882 Molecule *molecule=molecules.Current(ithread);
00883 int type=molecule->Type();
00884 if(type!=VOIDSPECIE&&molecule->index()%Run::nthreads==ithread) {
00885 REAL *v=molecule->Velocity();
00886 for(int i=0;i<DIM;i++) {
00887 if(v[i]>SMALL) {
00888 REAL dt=maxjump/v[i];
00889 if(dt<Run::time.step&&dt>SMALL)
00890 #pragma omp critical
00891 {
00892 Run::time.step=dt;
00893 }
00894 }
00895 }
00896 }
00897 molecules.goNext(ithread);
00898 } while(!molecules.isFirst(ithread));
00899 }
00900 Run::time.next=Run::time.current+Run::time.step;
00901 int nmolecules=0;
00902
00903 while(molecules.number()>0) {
00904 molecules.goFirst();
00905
00906 molecules.goFirst();
00907 do {
00908 Molecule *molecule=molecules.Current();
00909 int type=molecule->Type();
00910 if(type!=VOIDSPECIE&&molecule->Time()<Run::time.next) {
00911 #ifdef ADSORPTION
00912 if(molecule->ReleaseTime()<=currenttime) {
00913 #endif
00914 molecule->Move();
00915 int newmolecules=boundary(molecule);
00916 type=molecule->Type();
00917 if(type==VOIDSPECIE) {
00918 Container<Molecule> *gridmolecules=GridContainer::nodes+molecule->GridCell();
00919 Ptr<Molecule> *ptr=molecule->getPtr();
00920 if(!gridmolecules->remove(ptr)) {
00921 cout<<"Failed to remove node from the grid array\n";
00922 exit(1);
00923 }
00924 molecule->putPtr(NULL);
00925 molecules.remove0();
00926 if(molecules.number()==0) break;
00927 } else {
00928
00929 GridContainer::put(molecule);
00930
00931
00932
00933 molecules.goNext();
00934 }
00935 #ifdef ADSORPTION
00936 }
00937 #endif
00938 } else
00939 molecules.goNext();
00940 } while(!molecules.isFirst());
00941
00942 molecules.goFirst();
00943 int nmoved=0;
00944 nmolecules=0;
00945 do {
00946 Molecule *molecule=molecules.Current();
00947 int type=molecule->Type();
00948 if(type==VOIDSPECIE) {
00949 Container<Molecule> *gridmolecules=GridContainer::nodes+molecule->GridCell();
00950 Ptr<Molecule> *ptr=molecule->getPtr();
00951 if(!gridmolecules->remove(ptr)) {
00952 cout<<"Failed to remove node from the grid array\n";
00953 exit(1);
00954 }
00955 molecule->putPtr(NULL);
00956 molecules.remove0();
00957 if(molecules.number()==0)break;
00958 } else {
00959
00960 if(molecule->Time()>=Run::time.next)nmoved++;
00961 molecule->index(nmolecules++);
00962 molecules.goNext();
00963 }
00964 } while(!molecules.isFirst());
00965 if(nmoved>=nmolecules)break;
00966
00967
00968 #pragma omp parallel
00969 {
00970 int ithread = omp_get_thread_num();
00971 molecules.goFirst(ithread);
00972 do {
00973 Molecule *molecule=molecules.Current(ithread);
00974 int type=molecule->Type(),imolecule=molecule->index();
00975 if(type!=VOIDSPECIE&&imolecule%Run::nthreads==ithread) {
00976 REAL time=molecule->Time();
00977 if(time<Run::time.next) {
00978 molecule->TimeStep(Run::time.next-time+SMALL);
00979 molecule->ready();
00980 } else {
00981 molecule->TimeStep(SMALL);
00982 molecule->done();
00983 }
00984 }
00985 molecules.goNext(ithread);
00986 } while(!molecules.isFirst(ithread));
00987 interaction(ithread);
00988 }
00989 if(nmoved>=nmolecules)break;
00990 }
00991
00992 for(int i=0;i<nspecies;i++) {
00993 distribution[i]=0;
00994 }
00995 REAL
00996 energy=0.0,
00997 totalKineticEnergy=0.0,
00998 totalInternalEnergyPdof=0.0,
00999 totalInternalEnergy=0.0;
01000 molecules.goFirst();
01001 nmolecules=0;
01002 if(molecules.number()>0)
01003 do {
01004 Molecule *molecule=molecules.Current();
01005 int type=molecule->Type();
01006 if(type!=VOIDSPECIE) {
01007
01008 REAL
01009 cp=species[type].Cp(),
01010 dof=DOF(cp),
01011 dofi=DOFI(dof),
01012 internalEnergy=molecule->InternalEnergy();
01013 totalKineticEnergy+=molecule->KineticEnergy();
01014 totalInternalEnergy+=internalEnergy;
01015 totalInternalEnergyPdof+=dofi>numeric_limits<REAL>::epsilon()?internalEnergy/dofi:(REAL)0.0,
01016 distribution[type]++;
01017 nmolecules++;
01018 }
01019 molecules.goNext();
01020 } while(!molecules.isFirst());
01021 energy=totalInternalEnergy+totalKineticEnergy;
01022 REAL totalKineticEnergyPdof=onethird*totalKineticEnergy,
01023 totalEnergyPdof=totalKineticEnergyPdof+totalInternalEnergyPdof;
01024 if(Run::option.xterm) {
01025 const int ncols=6;
01026 printf(
01027 "%c[HTime=%-9.3e Step=%-9.3e Memory=%d%% N=%d T=%-6.2fK E=%-9.3eJ K=%-9.3eJ I=%-9.3eJ%c[K\n",
01028 ESC,currenttime,Run::time.step,
01029 (int)(round(100.0*(REAL)nmolecules/(REAL)molecules.maxnumber())),
01030 nmolecules,
01031 twothirds*AtomicMassUnit*totalKineticEnergy/(BoltzmannConstant*(REAL)nmolecules),
01032 AtomicMassUnit*energy,
01033 AtomicMassUnit*totalKineticEnergyPdof,
01034 AtomicMassUnit*totalInternalEnergyPdof,
01035 ESC
01036 );
01037 int ispecie=0,jspecie=0;
01038 while(ispecie<nspecies) {
01039 int icol=0;
01040 for(jspecie=ispecie;ispecie<nspecies;ispecie++){
01041 if(distribution[ispecie]>0) {
01042 if(icol++==ncols) break;
01043 printf("%%%-9s",species[ispecie].Id());
01044 }
01045 }
01046 if(ispecie<=nspecies)printf("%c[K\n",ESC);
01047 icol=0;
01048 for(;jspecie<ispecie;jspecie++){
01049 if(distribution[jspecie]>0) {
01050 if(icol++==ncols) break;
01051 printf("%-9.5f ",100.0*(REAL)distribution[jspecie]/(REAL)molecules.number());
01052 }
01053 }
01054 if(ispecie<=nspecies)printf("%c[K\n",ESC);
01055 }
01056 fflush(stdout);
01057 }
01058 currenttime=Run::time.next;
01059 Run::time.current=currenttime;
01060 injection();
01061 if(Run::time.output>0.0 && (int)floor(currenttime/Run::time.output)>IO::ioutput) {
01062 IO::ioutput=(int)floor(currenttime/Run::time.output);
01063 if(Run::option.verbose) cout<<"Writing output set "<<Run::outputname<<endl;
01064 save(Run::outputname);
01065 }
01066 }
01067 clock_t end_wtime=time((time_t) NULL);
01068 REAL elapsed_time=difftime(end_wtime,start_wtime);
01069 if (Run::option.debug)
01070 cout << " elapsed time: "<<elapsed_time<<"s, time per iteration: "<<elapsed_time/(REAL)niter<<"s\n";
01071 if(currenttime>endtime) currenttime=starttime+currenttime-endtime;
01072 if(currenttime<starttime) currenttime=endtime-starttime+currenttime;
01073 Run::time.current=currenttime;
01074 return iter;
01075 }
01076 int Domain::boundary(Molecule *a) {
01077 int newmolecules=0;
01078
01079 REAL dx=0.0,
01080
01081 *x=a->Coordinates(),
01082 *v=a->Velocity();
01083 int idir=-1, iside=-1;
01084 for(int i=0;i<DIM;i++) {
01085 if(x[i]<=XMIN[i]) {idir=i;iside=0;break;}
01086 if(x[i]>=XMAX[i]) {idir=i;iside=1;break;}
01087 }
01088 if(idir<0) return insideBoundary;
01089 dx=XMAX[idir]-XMIN[idir];
01090 Boundary *boundary=&boundaries[idir][iside];
01091 enum BoundaryTypes btype=boundary->type;
01092 switch(btype) {
01093 case periodicBoundary:
01094
01095 if(iside==0) x[idir]+=dx;
01096 else x[idir]-=dx;
01097
01098
01099
01100
01101
01102 break;
01103
01104 case elasticBoundary:
01105 {
01106 v[idir]=-v[idir];
01107 x[idir]=2.0*bounds[iside][idir]-x[idir];
01108
01109 int jdir=-1, jside=-1;
01110 for(int i=1;i<DIM;i++) {
01111 int j=(idir+i)%DIM;
01112 if(x[j]<=bounds[0][j]){jdir=j;jside=0;break;}
01113 if(x[j]>=bounds[1][j]){jdir=j;jside=1;break;}
01114 }
01115 if(jdir>=0){
01116 v[jdir]=-v[jdir];
01117 x[jdir]=2.0*bounds[jside][jdir]-x[jdir];
01118 }
01119 }
01120 break;
01121 case openBoundary:
01122 a->Type(VOIDSPECIE);
01123 return btype;
01124 default:
01125 cerr<<"Unknown BoundaryType "<<btype<<" encountered\n";
01126 }
01127
01128 int ispecie=a->Type();
01129 Reaction *reaction=boundary->reactions+ispecie;
01130 Reaction::Outcome *outcome=reaction->First();
01131 if(outcome==NULL) return btype;
01132 REAL r=RND,
01133 probability=0.0;
01134 bool reacted=false;
01135 do {
01136 probability+=outcome->Probability();
01137 if(r<probability) {
01138 reacted=true;
01139 REAL
01140 mAold=species[ispecie].Mass(),
01141 enthalpy=outcome->Enthalpy(),
01142 #ifdef ADSORPTION
01143 releasetime=outcome->Time(),
01144 #endif
01145 energy=a->KineticEnergy()+a->InternalEnergy()+enthalpy;
01146 int
01147 typeA=outcome->Product(0),
01148 typeB=outcome->Product(1);
01149 if(typeA<0||typeA>=VOIDSPECIE) {
01150 cerr<<"Can't locate specie type "<<typeA<<" for boundary reaction\n";
01151 exit(1);
01152 }
01153 Specie *sA=species+typeA;
01154 REAL
01155 mA=sA->Mass(),
01156 cpA=sA->Cp(),
01157
01158
01159
01160
01161
01162
01163
01164 dofA=DOF(cpA);
01165 REAL dofiA=DOFI(dofA),
01166 vA=LENGTH(v),
01167 keA=0.0,ieA=0.0;
01168 if(typeB!=VOIDSPECIE) {
01169 Specie *sB=species+typeB;
01170 REAL
01171 mB=sB->Mass(),
01172 cpB=sB->Cp(),
01173
01174
01175
01176
01177
01178
01179
01180 dofB=DOF(cpB);
01181 REAL dofiB=DOFI(dofB),
01182 ke=energy>0.0?2*DOFK/(dofA+dofB)*energy:0.0,
01183 ie=energy-ke;
01184 keA=0.5*ke;
01185 REAL
01186 keB=keA,
01187 dofi=dofiA+dofiB,
01188 ieB=0.0;
01189 if(dofi>SMALL) {
01190 ieA=dofiA/dofi*ie;
01191 ieB=ie-ieA;
01192 }
01193 REAL
01194 m=mA+mB,
01195 vAnew=sqrt(2.0*keA/mA),
01196 vBnew=sqrt(2.0*keB/mB);
01197 if(vA>SMALL) {
01198 REAL c=vAnew/vA;
01199 MULC(DIM,v,c);
01200 } else {
01201 v[idir]=iside==0?vAnew:-vAnew;
01202 }
01203
01204 Molecule *b=molecules.insertAfter();
01205
01206 REAL
01207 *u=b->Velocity(),
01208 *y=b->Coordinates();
01209 for(int i=0;i<DIM;i++){
01210 y[i]=x[i];
01211 u[i]=0.0;
01212 }
01213 REAL d=RND*vBnew*Run::time.step+0.5*(sA->Size()+sB->Size());
01214 if(iside==0) {
01215 u[idir]=vBnew;
01216 x[idir]+=sA->Size();
01217 y[idir]+=d;
01218 } else {
01219 u[idir]=-vBnew;
01220 x[idir]-=sA->Size();
01221 y[idir]-=d;
01222 }
01223 b->InternalEnergy(ieB);
01224 b->Type(typeB);
01225 b->Time(Run::time.current);
01226 b->TimeStep(Run::time.next-Run::time.current+SMALL);
01227 b->ready();
01228 #ifdef ADSORPTION
01229 if(residencetime>SMALL) q->ReleaseTime(Run::time.current+time);
01230 #endif
01231 newmolecules++;
01232 } else {
01233
01234 keA=energy>0.0?DOFK/dofA*energy:0.0;
01235 ieA=energy-keA;
01236 REAL vAnew=sqrt(2.0*keA/mA);
01237 REAL d=RND*vAnew*Run::time.step+(sA->Size());
01238 if(vA>SMALL) {
01239 REAL c=vAnew/vA;
01240 MULC(DIM,v,c);
01241 } else {
01242 v[idir]=iside==0?vAnew:-vAnew;
01243 if(iside==0) {
01244 x[idir]+=d;
01245 } else {
01246 x[idir]-=d;
01247 }
01248 }
01249 }
01250 a->InternalEnergy(ieA);
01251 a->Type(typeA);
01252 #ifdef ADSORPTION
01253 if(residencetime>SMALL) a->ReleaseTime(Run::time.current+time);
01254 #endif
01255 }
01256 } while(reacted==false&&(outcome=reaction->Next())!=NULL);
01257 return newmolecules;
01258 }
01260 void Domain::injection() {
01261 REAL dt=(REAL)Run::time.step*1e-9;
01262
01263 for(int idir=0;idir<DIM;idir++)
01264 for(int iside=0;iside<2;iside++) {
01265 Boundary *boundary=&boundaries[idir][iside];
01266 REAL area=1.0e-18*(REAL)boundary->area,
01267 temperature=boundary->temperature;
01268 List<Gas> *gases=&boundary->gases;
01269 if(gases->number()>0) {
01270 gases->goFirst();
01271 do {
01272 Gas *gas=gases->Current();
01273 Specie *specie=gas->specie;
01274
01279 REAL
01280 mu=(REAL)specie->Mass(),
01281 mass=AtomicMassUnit*mu,
01282 den=(REAL)gas->density,
01283 energy=BoltzmannConstant*temperature,
01284 Vavx=sqrt(energy/mass),
01285
01286 freq=den*Vavx/(2.0*mass),
01287 totinj=freq*area*dt,
01288 wholeinj=floor(totinj),
01289 fracinj=totinj-wholeinj;
01290
01291 Molecule *molecule=NULL;
01292 int ninj=(int)wholeinj;
01293 for(int inj=0;inj<ninj;inj++) {
01294 boundary->Inject((int)(specie-Species::species),Vavx,temperature,&molecules);
01295 }
01297 if(RND<=fracinj) {
01298 boundary->Inject((int)(specie-Species::species),Vavx,temperature,&molecules);
01299 }
01300 gases->goNext();
01301 } while(!gases->isFirst());
01302 }
01303 }
01304 }
01305 #ifdef LOCAL // ACCELERATION BY USING LOCAL INTERACTION MODEL: NO COLLISION DETECTION YET
01306 void Domain::interaction() {
01307 using namespace GridContainer;
01308
01309 int gridcell[DIM];
01310 REAL dx[DIM];
01311 for(int i=0;i<DIM;i++) dx[i]=XMAX[i]-XMIN[i];
01312 for(int ix=1;ix<ncells[0]-1;ix++) { gridcell[0]=ix;
01313 for(int iy=1;iy<ncells[1]-1;iy++) { gridcell[1]=iy;
01314 for(int iz=1;iz<ncells[2]-1;iz++) { gridcell[2]=iz;
01315 int icell=index(gridcell);
01316 Container<Molecule> *particles=nodes+icell;
01317 if(particles->number()==0) continue;
01318 particles->setFirstLast();
01319 do {
01320 Molecule *a=particles->First();
01321 if(a->Type()!=VOIDSPECIE) {
01322 Interaction status=missed;
01323
01324 if(particles->number()>1) {
01325 while (status==missed&&!particles->isLast()) {
01326 particles->goNext();
01327 Molecule *b=particles->Current();
01328 if(b->Type()!=VOIDSPECIE) {
01329 status=interact(a, b);
01330 if(status==annihalated) {
01331 if(b->Type()==VOIDSPECIE) {
01332 particles->remove(b->getPtr());
01333 b->putPtr(NULL);
01334 break;
01335 }
01336 else {
01337 cout<<"Molecule annihalated but not VOIDed!\n";cout.flush();exit(1);
01338 }
01339 }
01340 }
01341 }
01342 }
01343
01344 int neibcell[DIM];
01345 for(int jx=-1;status==missed&&jx<2;jx++) { neibcell[0]=ix+jx;
01346 for(int jy=-1;status==missed&&jy<2;jy++) { neibcell[1]=iy+jy;
01347 for(int jz=-1;status==missed&&jz<2;jz++) { neibcell[2]=iz+jz;
01348 if(jx==0&&jy==0&&jz==0) continue;
01349 int jcell=index(neibcell);
01350 Container<Molecule> *neibs=nodes+jcell;
01351 if(neibs->number()==0) continue;
01352 neibs->goFirst();
01353 do {
01354 Molecule *b=neibs->Current();
01355 if(b->Type()!=VOIDSPECIE) {
01356 status=interact(a, b);
01357 if(status==annihalated) {
01358 if(b->Type()==VOIDSPECIE) {
01359 neibs->remove(b->getPtr());
01360 b->putPtr(NULL);
01361
01362
01363 break;
01364 }
01365 else {
01366 cout<<"Molecule annihalated but not VOIDed!\n";cout.flush();exit(1);
01367 }
01368 }
01369 }
01370 neibs->goNext();
01371 } while(status==missed&&!neibs->isFirst());
01372 }
01373 }
01374 }
01375 }
01376 particles->FirstNext();
01377 } while(!particles->isFirstLast());
01378 particles->setFirst();
01379 }
01380 }
01381 }
01382 }
01383 void Domain::interaction(int icounter) {
01384 using namespace GridContainer;
01385
01386 int gridcell[DIM];
01387 REAL
01388 nexttime=Run::time.next,
01389 dx[DIM];
01390 for(int i=0;i<DIM;i++) dx[i]=XMAX[i]-XMIN[i];
01391 for(int ix=1;ix<ncells[0]-1;ix++) { gridcell[0]=ix;
01392 for(int iy=1;iy<ncells[1]-1;iy++) { gridcell[1]=iy;
01393 for(int iz=1;iz<ncells[2]-1;iz++) { gridcell[2]=iz;
01394 int icell=index(gridcell);
01395 Container<Molecule> *particles=nodes+icell;
01396 if(particles->number()==0) continue;
01397
01398 particles->setFirstLast(icounter);
01399 do {
01400 Molecule *a=particles->First(icounter);
01401 int iparticle=a->index();
01402 if(a->Type()!=VOIDSPECIE &&
01403 a->Time()<nexttime &&
01404 iparticle%Run::nthreads==icounter &&
01405 !a->isdone()) {
01406 a->done();
01407 Interaction status=missed;
01408 {
01409 Molecule *b=a->Target();
01410 if(b!=NULL&&!b->isdone()&&b->Type()!=VOIDSPECIE){
01411 if(!b->islocked()) {
01412
01413 b->lock();
01414 status=interact(a, b);
01415 if(status==annihalated&&b->Type()!=VOIDSPECIE)
01416 b->Type(VOIDSPECIE);
01417 b->unlock();
01418 if(status!=missed){
01419 particles->FirstNext(icounter);
01420
01421 continue;
01422 } else a->Target(NULL);
01423 }
01424 }
01425 }
01426
01427 if(particles->number()>1) {
01428 while (status==missed&&!particles->isLast(icounter)) {
01429 particles->goNext(icounter);
01430 Molecule *b=particles->Current(icounter);
01431 if(b->Type()!=VOIDSPECIE&&
01432 b->Time()<nexttime &&
01433 !b->isdone()) {
01434 if(!b->islocked()) {
01435
01436 b->lock();
01437 status=interact(a, b);
01438 if(status==annihalated&&b->Type()!=VOIDSPECIE)
01439 b->Type(VOIDSPECIE);
01440 b->unlock();
01441 }
01442 }
01443 }
01444 }
01445
01446 int neibcell[DIM];
01447 for(int jx=-1;status==missed&&jx<2;jx++) { neibcell[0]=ix+jx;
01448 for(int jy=-1;status==missed&&jy<2;jy++) { neibcell[1]=iy+jy;
01449 for(int jz=-1;status==missed&&jz<2;jz++) { neibcell[2]=iz+jz;
01450 if(jx==0&&jy==0&&jz==0) continue;
01451 int jcell=index(neibcell);
01452 Container<Molecule> *neibs=nodes+jcell;
01453 int nneibs=neibs->number();
01454 if(nneibs<=0) continue;
01455 neibs->goFirst(icounter);
01456 do {
01457 Molecule *b=neibs->Current(icounter);
01458 if(b->Type()!=VOIDSPECIE &&
01459 b->Time()<nexttime &&
01460 !b->isdone()) {
01461 if(!b->islocked()) {
01462
01463 b->lock();
01464 status=interact(a, b);
01465 if(status==annihalated&&b->Type()!=VOIDSPECIE)
01466 b->Type(VOIDSPECIE);
01467 b->unlock();
01468 }
01469 }
01470 neibs->goNext(icounter);
01471 nneibs--;
01472 } while(status==missed&&!neibs->isFirst(icounter)&&nneibs>0);
01473 }
01474 }
01475 }
01476 }
01477 particles->FirstNext(icounter);
01478
01479 } while(!particles->isFirstLast(icounter));
01480 particles->setFirst(icounter);
01481 }
01482 }
01483 }
01484 }
01485 #ifdef HARDBALLS
01486 Interaction Domain::interact(Molecule *a, Molecule *b) {
01487 Interaction status=missed;
01488 int typeA=a->Type(), typeB=b->Type();
01489 if(typeA==VOIDSPECIE||typeB==VOIDSPECIE){
01490 #ifdef DEBUG
01491 cerr<<"Wrong reactants species: \n"
01492 <<"\tA: type="<<typeA<<", id="<<species[typeA].Id()<<endl
01493 <<"\tB: type="<<typeB<<", id="<<species[typeB].Id()<<endl;
01494 #endif
01495 return annihalated;
01496 }
01497 Specie
01498 *spa=species+typeA,
01499 *spb=species+typeB;
01500 REAL
01501 *xa=a->Coordinates(),
01502 *va=a->Velocity(),
01503 *xb=b->Coordinates(),
01504 *vb=b->Velocity();
01505 REAL
01506 ma=spa->Mass(),
01507 da=spa->Size(),
01508 mb=spb->Mass(),
01509 db=spb->Size(),
01510 d0=0.5*(da+db),
01511 ta=a->Time(),
01512 tb=b->Time(),
01513 dtab=b->Time()-a->Time(),
01514 dvab2=0.0,
01515 dabt2=0.0,
01516 Dab[DIM],dab2=0.0;
01517 for(int i=0;i<DIM;i++)
01518 { REAL r=xb[i]-xa[i];
01519 Dab[i]=r;
01520 dab2+=r*r;
01521 }
01522
01523 for(int i=0;i<DIM;i++){
01524 REAL dv=vb[i]-va[i];
01525 dvab2+=dv*dv;
01526 }
01527 dabt2=dab2+dvab2*dtab*dtab;
01528 if(sqrt(dabt2)<d0) {
01529
01530 REAL
01531 m=ma+mb,mi=1.0/m,
01532 dab=sqrt(dab2), dabi=1.0/dab,
01533 d1[DIM],
01534 w[DIM],
01535 ua[DIM],ub[DIM],
01536
01537 una,unb;
01538
01539
01540
01541 for(int k=0;k<DIM;k++) {
01542 REAL
01543 vak=(REAL)va[k],
01544 vbk=(REAL)vb[k],
01545 wk=mi*(ma*vak+mb*vbk);
01546 ua[k]=vak-wk;
01547 ub[k]=vbk-wk;
01548 w[k]=wk;
01549 d1[k]=Dab[k]*dabi;
01550 }
01551
01552 una=SCLP(ua,d1);
01553 unb=SCLP(ub,d1);
01554
01555
01556
01557
01558
01559
01560 if(una-unb>SMALL) {
01561 REAL
01562 ua2old=SCLP(ua,ua),uaold=sqrt(ua2old),
01563 ub2old=SCLP(ub,ub),ubold=sqrt(ub2old),
01564 keold=0.5*(ma*ua2old+mb*ub2old),
01565 ke=keold,
01566 ieold=(REAL)a->InternalEnergy()+(REAL)b->InternalEnergy(),
01567 ie=ieold,
01568 energy=keold+ieold;
01569 status=collided;
01570
01571 for(int k=0;k<DIM;k++) {
01572 ua[k]-=2.0*d1[k]*una;
01573 ub[k]-=2.0*d1[k]*unb;
01574 }
01575
01576 Reaction *reaction=reactions+typeA*nspecies+typeB;
01577 Reaction::Outcome *outcome=reaction->First();
01578 #ifdef DEBUG
01579 int typeAold=typeA, typeBold=typeB;
01580 #endif
01581 if(outcome!=NULL) {
01582 REAL
01583 r=RND,
01584 probability=0.0;
01585 do {
01586 probability+=outcome->Probability();
01587 if( r < probability &&
01588 energy >= outcome->ActivationEnergy()
01589 ) {
01590 typeA=outcome->Product(0);
01591 typeB=outcome->Product(1);
01592 #ifdef DEBUG
01593 if(typeA==VOIDSPECIE) {
01594 cout<<"typeA="<<typeA<<endl;
01595 typeA=typeB; typeB=VOIDSPECIE;
01596 }
01597 #endif
01598 if(typeB!=VOIDSPECIE) {
01599 status=reacted;
01600 } else {
01601 status=annihalated;
01602 }
01603 a->Type(typeA);
01604 spa=species+typeA;
01605 b->Type(typeB);
01606 spb=species+typeB;
01607 break;
01608 }
01609 } while((outcome=reaction->Next())!=NULL);
01610 }
01611 REAL
01612 cpa=spa->Cp(),
01613 dofa=DOF(cpa),
01614 dofia=DOFI(dofa),
01615 cpb=spb->Cp(),
01616 dofb=DOF(cpb),
01617 dofib=DOFI(dofb),
01618 dof=dofa+dofb,
01619 dofk=2.0*DOFK,
01620 dofi=dofia+dofib,
01621 rdofia=0.0, rdofib=0.0;
01622 if(dofi>SMALL) {
01623 rdofia=dofia/dofi;
01624 rdofib=1.0-rdofia;
01625 }
01626 REAL
01627 dofkc=3.0,
01628 dofic=dofi,
01629 dofc=dofkc+dofic,
01630
01631 resa=1.0, resb=1.0;
01632 if(status>=reacted) {
01633 REAL
01634 manew=spa->Mass(),
01635 mbnew=spb->Mass(),
01636 enthalpy=outcome->Enthalpy();
01637 energy+=enthalpy;
01638 #ifdef DEBUG
01639 if(fabs(ma+mb-manew-mbnew)>=SMALL) {
01640 cerr<<"Mass conservation violation in interaction:\n"
01641 <<species[typeAold].Id()<<'('<<species[typeAold].Mass()<<") + "
01642 <<species[typeBold].Id()<<'('<<species[typeBold].Mass()<<") = "
01643 <<ma+mb<<" != "
01644 <<spa->Id()<<'('<<spa->Mass()<<") + "
01645 <<spb->Id()<<'('<<spb->Mass()<<") = "
01646 <<manew+mbnew<<endl;
01647 exit(1);
01648 }
01649 #endif
01650 if(status!=annihalated) {
01651 if(energy>0.0) {
01652
01653 ke=dofkc*energy/dofc;
01654
01655 ie=energy-ke;
01656 if(uaold>SMALL&&ubold>SMALL) {
01657 REAL
01658
01659
01660
01661 uanew=sqrt(ke/ma),
01662 ubnew=sqrt(ke/mb),
01663 uaratio=uanew/uaold,
01664 ubratio=ubnew/ubold;
01665 for(int i=0;i<DIM;i++) {
01666 ua[i]*=uaratio;
01667 ub[i]*=ubratio;
01668 }
01669 }
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682 a->InternalEnergy(rdofia*ie);
01683 b->InternalEnergy(rdofib*ie);
01684 } else {
01685 REAL
01686 ea=rdofia*enthalpy,
01687 eb=enthalpy-ea;
01688 a->InternalEnergy(a->InternalEnergy()+ea);
01689 b->InternalEnergy(b->InternalEnergy()+eb);
01690 }
01691 } else {
01692 for(int i=0;i<DIM;i++) { ua[i]=ub[i]=0.0; }
01693 a->InternalEnergy(energy);
01694 b->InternalEnergy(0.0);
01695 }
01696 } else {
01697
01698 ke=energy>0.0?dofkc*energy/dofc:0.0;
01699
01700 ie=energy-ke;
01701 if(uaold>SMALL&&ubold>SMALL) {
01702 REAL
01703
01704
01705
01706 uanew=sqrt(ke/ma),
01707 ubnew=sqrt(ke/mb),
01708 uaratio=uanew/uaold,
01709 ubratio=ubnew/ubold;
01710 for(int i=0;i<DIM;i++) {
01711 ua[i]*=uaratio;
01712 ub[i]*=ubratio;
01713 }
01714 }
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729 a->InternalEnergy(rdofia*ie);
01730 b->InternalEnergy(rdofib*ie);
01731
01732 }
01733
01734
01735
01736
01737
01738
01739 for(int k=0;k<DIM;k++) {
01740 va[k]=(REAL)(ua[k]+w[k]);
01741 vb[k]=(REAL)(ub[k]+w[k]);
01742 }
01743 }
01744
01745 a->Target(NULL);
01746 b->Target(NULL);
01747 }
01748 else {
01749 REAL
01750 *vb=b->Velocity(),va2=0.0,vb2=0.0,
01751 dv2=0.0,dvi,
01752 dV[DIM],
01753 Xa[DIM],Xb[DIM];
01754
01755 for(int i=0;i<DIM;i++){
01756 REAL
01757 va0=va[i],
01758 vb0=vb[i],
01759 dv=va0-vb0;
01760 dv2+=dv*dv;
01761 dV[i]=dv;
01762 Xa[i]=xa[i];
01763 Xb[i]=xb[i];
01764 va2+=va0*va0;
01765 vb2+=vb0*vb0;
01766 }
01767
01768 if(dtab>0.0)for(int i=0;i<DIM;i++)Xa[i]+=va[i]*dtab;
01769 else for(int i=0;i<DIM;i++)Xb[i]-=vb[i]*dtab;
01770 dab2=0.0;
01771 for(int i=0;i<DIM;i++){
01772 REAL r=Xb[i]-Xa[i];
01773 Dab[i]=r;
01774 dab2+=r*r;
01775 }
01776 REAL dab=sqrt(dab2);
01777 if(dv2>SMALL) {
01778 REAL
01779 d2=0.0,
01780 du=0.0,
01781 dx2=0.0,
01782 dvx=0.0,
01783 U[DIM];
01784 dvi=1.0/sqrt(dv2);
01785 for(int i=0;i<DIM;i++) {
01786 REAL u=dV[i]*dvi;
01787 du+=Dab[i]*u;
01788 U[i]=u;
01789 }
01790 for(int i=0;i<DIM;i++) {
01791 REAL
01792 r=Dab[i],
01793 dx=r-SMALL,
01794 dv=dV[i],
01795 du=U[i]*dab,
01796 di=r+du;
01797 d2+=di*di;
01798 dx2=+dx*dx;
01799 dvx=+dx*dv;
01800 }
01801 if(d0-SMALL>sqrt(d2)&&dvx>SMALL) {
01802 REAL dt=dx2/dvx,
01803 dta=a->TimeStep(),
01804 tb=b->Time(),
01805 tcolb=tb+b->TimeStep(),
01806 dtamin=0.5*species[typeA].Size()/sqrt(va2),
01807 dtbmin=0.5*species[typeB].Size()/sqrt(vb2);
01808 if(dtab>0.0)dt+=dtab;
01809
01810 REAL coltime=a->Time()+dt,
01811 dtb=coltime-tb;
01812 if(dt<a->TimeStep()&&dt>dtamin&&coltime<tcolb&&dtb>=0.0&&dt>dtbmin)
01813 #pragma omp critical //DDD should declare a protected instead
01814 {
01815
01816 a->TimeStep(dt);
01817 b->TimeStep(dtb);
01818 a->Target(b);
01819 b->Target(a);
01820 }
01821 }
01822 }
01823 }
01824 return status;
01825 }
01826 #else //NOT-HARD-BALLS: Lennard-Jones:
01827 bool Domain::interact(Molecule *a, Molecule *b) {
01828 bool collision=false;
01829 REAL
01830
01831 *xa=a->Coordinates(),
01832
01833 *xb=b->Coordinates(),
01834 cutoffdistance=Potential::Cutoff(),
01835 d[DIM],d0,d2=0.0;
01836 for(int i=0;i<DIM;i++)
01837 { REAL r=xb[i]-xa[i];
01838 d[i]=r;
01839 d2+=r*r;
01840 }
01841 d0=sqrt(d2);
01842 if(d0<cutoffdistance) {
01843
01844
01845
01846
01847 REAL
01848 dt=Run::time.step,
01849 f=Potential::force(d0),
01850 ma=species[a->Type()].Mass(),
01851 mb=species[b->Type()].Mass(),
01852 dva=f/ma*dt,
01853 dvb=f/mb*dt,
01854 *va=a->Velocity(),
01855 *vb=b->Velocity();
01856 for(int k=0;k<DIM;k++) {
01857 REAL r=d[k]/d0;
01858 va[k]+=dva*r;
01859 vb[k]+=dvb*r;
01860 }
01861 collision=true;
01862 }
01863 return collision;
01864 }
01865
01866 #endif
01867 #else // SIMPLE N^2 INTERACTION ALGORITHM:
01868 void Domain::interaction() {
01869
01870 REAL dx[DIM];
01871 for(int i=0;i<DIM;i++) dx[i]=XMAX[i]-XMIN[i];
01872 molecules.setFirstLast();
01873 do {
01874 Molecule *a=molecules.First();
01875 REAL
01876 ma=a->Mass(),
01877 ra=0.5*a->Size(),
01878 *va=a->Velocity(),
01879 xa[2][DIM];
01880 for(int i=0;i<DIM;i++) {
01881 REAL x=a->Coordinate(i);
01882 xa[0][i]=x;
01883 if(x-XMIN[i]<0.5*dx[i])
01884 xa[1][i]=x+dx[i];
01885 else
01886 xa[1][i]=x-dx[i];
01887 }
01888 molecules.goFirstNext();
01889 while (!molecules.isLast()) {
01890 Molecule *b=molecules.Current();
01891 REAL
01892 mb=b->Mass(),
01893 rb=0.5*b->Size(),
01894 *xb=b->Coordinates();
01895 for(int j=0;j<2;j++) {
01896 REAL d[DIM],d0,d2=0.0;
01897 for(int i=0;i<DIM;i++)
01898 { REAL r=xb[i]-xa[j][i];
01899 d[i]=r;
01900 d2+=r*r;
01901 }
01902 d0=sqrt(d2);
01903 if(d0<ra+rb) {
01904 REAL m=ma+mb,
01905 *vb=b->Velocity(),
01906 d1[DIM],
01907 w[DIM],
01908 ua[DIM],ub[DIM],
01909
01910 una,unb;
01911
01912 for(int k=0;k<DIM;k++)
01913 { w[k]=(ma*va[k]+mb*vb[k])/m;
01914 ua[k]=va[k]-w[k];
01915 ub[k]=vb[k]-w[k];
01916 d1[k]=d[k]/d0;
01917 }
01918 una=SCLP(ua,d1);
01919 unb=SCLP(ub,d1);
01920 if(una-unb>0.0) {
01921 for(int k=0;k<DIM;k++) {
01922 va[k]=ua[k]-2.0*d1[k]*una+w[k];
01923 vb[k]=ub[k]-2.0*d1[k]*unb+w[k];
01924 }
01925 }
01926 }
01927 }
01928 molecules.goNext();
01929 }
01930 molecules.FirstNext();
01931 } while(!molecules.isFirstLast());
01932 molecules.setFirst();
01933 }
01934 #endif
01935 void Domain::save(char *name) {
01936 char filename[MAXLINLEN];
01937 int nmolecules=molecules.number();
01938 if(nmolecules<=0) {
01939 if(Run::option.verbose) cout<< "No molecules saved\n";
01940 return;
01941 }
01942 sprintf(filename,"%s-%d.dat.gz",name,IO::ioutput);
01943 cout << "Saving task '"<<name<<"' into file: " << filename <<". Elapsed world time: "<<Run::elapsed()<<" sec"<<endl;
01944 gzFile out;
01945 if((out = gzopen(filename, "wb"))==NULL)
01946 { cerr << "CAN'T OPEN "<<filename<<endl;cerr.flush();
01947 exit(1);
01948 }
01949 if(Run::option.verbose||Run::option.debug)
01950 cout<<"Writing "<<nmolecules<<" molecules to "<<filename<<" ... ";
01951 gzprintf(out,"# Data output %d at time %lg for number of molecules %d\n",IO::ioutput,Run::time.current,nmolecules);
01952 if(Run::option.verbose||Run::option.debug)cout<<"Number of molecules: "<<nmolecules<<endl;
01953 int imolecule=0;
01954 molecules.goFirst();
01955 do {
01956 Molecule *molecule=molecules.Current();
01957 if(molecule->Type()!=VOIDSPECIE) {
01958 gzprintf(out,"%s",species[molecule->Type()].Id());
01959 for(int i=0;i<DIM;i++) gzprintf(out,"\t%lg",molecule->Coordinate(i));
01960 for(int i=0;i<DIM;i++) gzprintf(out,"\t%lg",molecule->Velocity(i));
01961 gzprintf(out,"\t%lg\t%lg\t%lg\n",molecule->InternalEnergy(),molecule->Time(),molecule->TimeStep());
01962
01963 }
01964 molecules.goNext();
01965 imolecule++;
01966 } while(!molecules.isFirst());
01967 gzprintf(out,"\n");
01968 if(imolecule!=nmolecules) {
01969 cerr<<"Number of molecules specified ("<<nmolecules
01970 <<") does not agree with actually found ("<<imolecule<<")\n";
01971 cerr.flush();exit(1);
01972 }
01973 if(Run::option.verbose) cout << molecules.number() << " molecules written\n";
01974 gzclose(out);
01975
01976 }
