ComboCode.cc.modeling.codes.Gastronoom

25 26 """ 27 Class that includes all methods required for creating a GATRoNOoM model. 28 29 """ 30

31 - def __init__(self,path_gastronoom='runTest',vic=None,sphinx=0,\ 32 replace_db_entry=0,cool_db=None,ml_db=None,sph_db=None,\ 33 skip_cooling=0,recover_sphinxfiles=0,\ 34 new_entries=[],single_session=0):

35 36 """ 37 Initializing an instance of a GASTRoNOoM modeling session. 38 39 A single cooling model and a given set of molecules and transitions are 40 calculated here or are retrieved from the databases. 41 42 @keyword path_gastronoom: modeling folder in GASTRoNOoM home 43 44 (default: 'runTest') 45 @type path_gastronoom: string 46 @keyword vic: the vic manager for running sphinx models on VIC3 47 48 (default: None) 49 @type vic: Vic() 50 @keyword sphinx: Running Sphinx? 51 52 (default: 0) 53 @type sphinx: bool 54 @keyword replace_db_entry: replace an entry in the databases with a 55 newly calculated model with a new model id 56 (e.g. if some general data not included in 57 the inputfiles is changed) 58 59 (default: 0) 60 @type replace_db_entry: bool 61 @keyword new_entries: The new model_ids when replace_db_entry is 1 62 of other models in the grid. These are not 63 replaced! 64 65 (default: []) 66 @type new_entries: list[str] 67 @keyword skip_cooling: Skip running cooling in case a model is not 68 found in the database, for instance if it is 69 already known that the model will fail 70 71 (default: 0) 72 @type skip_cooling: bool 73 @keyword recover_sphinxfiles: Try to recover sphinx files from the disk 74 in case they were correctly calculated, 75 but not saved to the database for one 76 reason or another. 77 78 (default: 0) 79 @type recover_sphinxfiles: bool 80 @keyword cool_db: the cooling database 81 82 (default: None) 83 @type cool_db: Database() 84 @keyword ml_db: the mline database 85 86 (default: None) 87 @type ml_db: Database() 88 @keyword sph_db: the sphinx database 89 90 (default: None) 91 @type sph_db: Database() 92 @keyword single_session: If this is the only CC session. Speeds up db 93 check. 94 95 (default: 0) 96 @type single_session: bool 97 98 """ 99 100 super(Gastronoom,self).__init__(code='GASTRoNOoM',\ 101 path=path_gastronoom,\ 102 replace_db_entry=replace_db_entry,\ 103 new_entries=new_entries,\ 104 single_session=single_session) 105 #-- Convenience path 106 cc.path.gout = os.path.join(cc.path.gastronoom,self.path) 107 self.vic = vic 108 self.sphinx = sphinx 109 cool_keys = os.path.join(cc.path.aux,'Input_Keywords_Cooling.dat') 110 ml_keys = os.path.join(cc.path.aux,'Input_Keywords_Mline.dat') 111 sph_keys = os.path.join(cc.path.aux,'Input_Keywords_Sphinx.dat') 112 self.cooling_keywords = [line.strip() 113 for line in DataIO.readFile(cool_keys) 114 if line] 115 self.mline_keywords = [line.strip() 116 for line in DataIO.readFile(ml_keys) 117 if line] 118 self.sphinx_keywords = [line.strip() 119 for line in DataIO.readFile(sph_keys) 120 if line] 121 DataIO.testFolderExistence(os.path.join(cc.path.gout,'data_for_mcmax')) 122 self.trans_bools = [] 123 self.mline_done = False 124 self.cool_done = False 125 #-- If a cooling model is in progress, the model manager will hold until 126 # the other cc session is finished. 127 self.in_progress = False 128 #-- Same for molecules. 129 self.molec_in_progress = [] 130 #-- Transitions in progress are simply ignored. 131 self.trans_in_progress = [] 132 self.cooling_molec_keys = ['ENHANCE_ABUNDANCE_FACTOR',\ 133 'ABUNDANCE_FILENAME',\ 134 'NUMBER_INPUT_ABUNDANCE_VALUES',\ 135 'KEYWORD_TABLE','MOLECULE_TABLE',\ 136 'ISOTOPE_TABLE'] 137 self.no_ab_molecs = ['12C16O','13C16O','1H1H16O','p1H1H16O',\ 138 '1H1H17O','p1H1H17O','1H1H18O','p1H1H18O'] 139 140 #- Read standard input file with all parameters that should be included 141 filename = os.path.join(cc.path.aux,'inputGASTRoNOoM.dat') 142 self.standard_inputfile = DataIO.readDict(filename,\ 143 comment_chars=['#','!']) 144 self.skip_cooling = skip_cooling 145 self.recover_sphinxfiles = recover_sphinxfiles 146 self.cool_db = cool_db 147 self.ml_db = ml_db 148 self.sph_db = sph_db

149 #self.pacs_db = pacs_db 150 151 152

153 - def addTransInProgress(self,trans):

154 155 ''' 156 Remember a transition in progress. They will be checked at the 157 end of a VIC run to see if they have been correctly calculated. If being 158 calculated in another CC session, they will not be checked. But the user 159 can always reload. 160 161 @param trans: The transition in progress 162 @type trans: Transition() 163 164 ''' 165 166 self.trans_in_progress.append(trans)

167 168 169

170 - def addMolecInProgress(self,molec):

171 172 ''' 173 Remember a molecule in progress. The model manager will wait until it 174 is finished. 175 176 @param molec: The molecule in progress 177 @type molec: Molecule() 178 179 ''' 180 181 self.molec_in_progress.append(molec)

182 183 184

185 - def execGastronoom(self,filename,subcode):

186 187 ''' 188 Execute a subprocess of GASTRoNOoM: cooling, mline or sphinx 189 190 @param filename: the full path+filename of the inputfile 191 @type filename: string 192 @param subcode: one of ['cooling','mline','sphinx'], the code to run 193 @type subcode: string 194 195 ''' 196 197 print '** Running %s...'%subcode 198 if not subcode.lower() in ['cooling','mline','sphinx']: 199 raise IOError('Subcode of GASTRoNOoM wrongly specified.') 200 subprocess.call(['echo %s | %s'%(filename,subcode.lower())],shell=True) 201 print '** DONE!' 202 print '***********************************'

203 204 205

206 - def makeIdLog(self,new_id,molec_id=None):

207 208 ''' 209 Make a text file with the original cooling id in it. 210 211 This is used when creating a Transition() from a sphinx outputfilename. 212 213 @param new_id: the new id for either mline or sphinx 214 @type new_id: string 215 @keyword molec_id: if given, an mline_id.log will be made including the 216 mline id to which molecules are linked 217 218 (default: None) 219 @type molec_id: string 220 221 ''' 222 223 DataIO.writeFile(filename=os.path.join(cc.path.gout,'models',new_id,\ 224 'cooling_id.log'),input_lines=[self.model_id]) 225 if not molec_id is None: 226 DataIO.writeFile(filename=os.path.join(cc.path.gout,'models',\ 227 new_id,'mline_id.log'),input_lines=[molec_id])

228 229 230

231 - def updateModel(self,model_id=None):

232 233 """ 234 Updating model name in command list, eg in case mline or sphinx require 235 a new model_id, different from the cooling model_id. 236 237 @keyword model_id: if None the instance's model_id is taken, otherwise 238 this one is used 239 240 (default: None) 241 @type model_id: string 242 243 """ 244 245 for par in ['OUTPUT_DIRECTORY','PARAMETER_FILE','OUTPUT_SUFFIX']: 246 self.command_list[par] = self.command_list[par].\ 247 replace(self.command_list[par][self.command_list[par].\ 248 find('model_'):self.command_list[par].find('model_')+25],\ 249 not model_id is None and model_id or self.model_id,2)

250 251 252

253 - def deleteCoolingId(self,model_id):

254 255 ''' 256 Delete everything from all databases that has anything to do with 257 this cooling id. 258 259 The databases on the hard disk are updated as well! 260 261 @param model_id: the cooling model_id 262 @type model_id: string 263 264 ''' 265 266 print 'Replacing database entries for old cooling id %s.'%model_id 267 del self.cool_db[model_id] 268 try: 269 del self.ml_db[model_id] 270 except KeyError: 271 pass 272 try: 273 del self.sph_db[model_id] 274 except KeyError: 275 pass

276 277 278

279 - def checkCoolingDatabase(self,molec_dict):

280 281 """ 282 Checking cooling database. 283 284 @param molec_dict: molecule info for this cooling model, ie CO and H2O 285 @type molec_dict: dict() 286 287 @return: The presence of the cooling model in the database 288 @rtype: bool 289 290 """ 291 292 #-- Lock the cooling database by opening it in read mode. It's closed 293 # once the database check is finalised. Note that in a case of a crash 294 # during the for loop, the python shell must be exited to unlock the 295 # sphinx database again. The sync() is now done only once at the very 296 # end since the file on the disk will not change. 297 if not self.single_session: self.cool_db.sync() 298 cool_dbfile = self.cool_db._open('r') 299 for i,model_id in enumerate(sorted(self.cool_db.keys())): 300 cool_dict = self.cool_db[model_id] 301 model_bool = self.cCL(self.command_list.copy(),cool_dict,'cooling',\ 302 extra_dict=molec_dict) 303 if model_bool: 304 if cool_dict.has_key('IN_PROGRESS'): 305 self.in_progress = True 306 print 'Cooling model is currently being calculated in a ' +\ 307 'different CC modeling session with ID %s.'\ 308 %(model_id) 309 self.model_id = model_id 310 self.updateModel() 311 finished = 1 312 break 313 elif self.replace_db_entry and model_id not in self.new_entries: 314 self.deleteCoolingId(model_id) 315 print 'Deleting cooling model with ID %s from database.'\ 316 %(model_id) + 'Calculating anew with id %s.'\ 317 %self.model_id 318 finished = 0 319 break 320 else: 321 print 'GASTRoNOoM cooling model has been calculated ' + \ 322 'before with ID %s.'%model_id 323 self.model_id = model_id 324 self.updateModel() 325 finished = 1 326 break 327 328 #-- Reached the end of db without match. Make new entry in db, in 329 # progress. Cant combine this with next line in case the last 330 # model gives a match. 331 if i == len(self.cool_db)-1: 332 print 'No match found in GASTRoNOoM cooling database. ' + \ 333 'Calculating new model.' 334 finished = 0 335 336 #-- In case of an empty db, the above loop is not accessed. 337 if not self.cool_db.keys(): 338 print 'No match found in GASTRoNOoM cooling database. ' + \ 339 'Calculating new model.' 340 finished = 0 341 342 #-- Add the model in progress to the cooling db 343 if finished == 0: 344 #-- OK to update. molec_dict is a copy. Because have to 345 # add the other input keywords for cooling to the H2O info. 346 # This is saved to the db 347 molec_dict.update(self.command_list) 348 molec_dict['IN_PROGRESS'] = 1 349 self.cool_db[self.model_id] = molec_dict 350 351 #-- Synchronize and unlock db. 352 cool_dbfile.close() 353 if not self.single_session: self.cool_db.sync() 354 return finished

355 356 357

358 - def checkMlineDatabase(self):

359 360 """ 361 Check mline database. 362 363 @return: The presence of mline models in the database, equivalent with 364 self.molec_list 365 @rtype: list[bool] 366 367 """ 368 369 #-- Lock the mline database by opening it in read mode. It's closed 370 # once the database check is finalised. Note that in a case of a crash 371 # during the for loop, the python shell must be exited to unlock the 372 # sphinx database again. The sync() is now done only once at the very 373 # end since the file on the disk will not change. 374 if not self.single_session: self.ml_db.sync() 375 ml_dbfile = self.ml_db._open('r') 376 if not self.ml_db.has_key(self.model_id): 377 self.ml_db[self.model_id] = dict([(self.model_id,dict())]) 378 for molec in self.molec_list: 379 molec.setModelId(self.model_id) 380 #-- Add an in-progress entry to the db 381 md = molec.makeDict(in_progress=1) 382 self.ml_db[self.model_id][self.model_id][molec.molecule] = md 383 self.ml_db.addChangedKey(self.model_id) 384 ml_dbfile.close() 385 self.ml_db.sync() 386 return [False]*len(self.molec_list) 387 388 model_bools = [] 389 for molec in self.molec_list: 390 for molec_id in [k for k,v in sorted(self.ml_db[self.model_id].items()) 391 if molec.molecule in v.keys()]: 392 db_molec_dict = self.ml_db[self.model_id][molec_id]\ 393 [molec.molecule] 394 if self.cCL(this_list=molec.makeDict(),\ 395 modellist=db_molec_dict,\ 396 code='mline',\ 397 ignoreAbun=molec.molecule in self.no_ab_molecs): 398 molec.setModelId(molec_id) 399 model_bools.append(True) 400 if db_molec_dict.has_key('IN_PROGRESS'): 401 self.addMolecInProgress(molec) 402 print 'Mline model is currently being ' + \ 403 'calculated in a different CC modeling '+\ 404 'session for %s with ID %s.'\ 405 %(molec.molecule,molec.getModelId()) 406 else: 407 print 'Mline model has been calculated before for '\ 408 '%s with ID %s.'\ 409 %(molec.molecule,molec.getModelId()) 410 break 411 412 #-- No match found. Calculate molecule anew. Now give it an id 413 # matching existing ids in terms of pars, or make a new one 414 if molec.getModelId() is None: 415 model_bools.append(False) 416 cool_dbd = self.ml_db[self.model_id] 417 for i,(k,v) in enumerate(sorted(cool_dbd.items())): 418 #-- No molecule in this id yet. 419 # This should NEVER occur! Raise error. 420 # Such empty dicts are removed upon fail/success check 421 # in doMline. 422 if not v.keys(): 423 raise KeyError('Empty molec id found in the database.'+\ 424 ' This should not be possible.') 425 426 #-- molecule already in this id. Keep searching, unless it 427 # is the last entry. Then we have to make a new id 428 if molec.molecule in v.keys() \ 429 and i != len(cool_dbd.keys())-1: 430 continue 431 432 #-- Now check if par dict is same except the molecule 433 # definition and some molecule specific parameters. Only 434 # check this specific keys list. If match, choose this ID 435 # Make sure molec is not in v.keys() already (to avoid 436 # the case of this being the last entry of the db) 437 if not molec.molecule in v.keys(): 438 mm = v.keys()[0] 439 cks = ['OUTER_R_MODE','CHANGE_DUST_TO_GAS_FOR_ML_SP',\ 440 'DUST_TO_GAS_CHANGE_ML_SP','STARFILE',\ 441 'USE_STARFILE','USE_NO_MASER_OPTION',\ 442 'USE_MASER_IN_SPHINX','FEHLER','N_FREQ',\ 443 'START_APPROX','USE_FRACTION_LEVEL_CORR',\ 444 'FRACTION_LEVEL_CORR','NUMBER_LEVEL_MAX_CORR'] 445 if self.cCL(this_list=molec.makeDict(),\ 446 modellist=v[mm],code='mline',\ 447 check_keys=cks): 448 break 449 450 #-- If end of dict is reached, no proper match was found 451 # Create new id and put it in k. 452 if i == len(cool_dbd.keys())-1: 453 #-- First check if the original model id is already 454 # in use (it may have been deleted if first use 455 # failed). Never copy output for this though, 456 # since it is the original 457 if not cool_dbd.has_key(self.model_id): 458 k = self.model_id 459 else: 460 k = self.makeNewId() 461 self.makeIdLog(new_id=k) 462 self.copyOutput(molec,self.model_id,k) 463 self.ml_db[self.model_id][k] = dict() 464 465 #-- It is possible cool_dbd is actually empty. Then use the 466 # model id as molec id. No need for an id log in this case 467 # nor any copying of output. 468 if not cool_dbd.keys(): 469 k = self.model_id 470 self.ml_db[self.model_id][self.model_id] = dict() 471 472 #-- The last k value (or the new one in case end of dict was 473 # reached) is set as the model id. 474 molec.setModelId(k) 475 476 #-- Add an in-progress entry to the db 477 md = molec.makeDict(in_progress=1) 478 self.ml_db[self.model_id][k][molec.molecule] = md 479 self.ml_db.addChangedKey(self.model_id) 480 481 #-- Inform the user 482 print 'Mline model for %s '%molec.molecule + \ 483 'has not been calculated before. Calculate anew with '+ \ 484 'ID %s.'%(molec.getModelId()) 485 ml_dbfile.close() 486 if not self.single_session: self.ml_db.sync() 487 return model_bools

488 489 490

491 - def checkSphinxDatabase(self):

492 493 """ 494 Check Sphinx database. 495 496 The presence of the sphinx models in the database is saved in 497 self.trans_bools, equivalent to self.trans_list 498 499 """ 500 501 #-- Remember which molecules have been added to new id, if applicable 502 copied_molecs = [] 503 504 #-- Lock the sphinx database by opening it in read mode. It's closed 505 # once the database check is finalised. Note that in a case of a crash 506 # during the for loop, the python shell must be exited to unlock the 507 # sphinx database again. The sync() is now done only once at the very 508 # end since the file on the disk will not change. 509 if not self.single_session: self.sph_db.sync() 510 sph_dbfile = self.sph_db._open('r') 511 if not self.sph_db.has_key(self.model_id): 512 self.sph_db[self.model_id] = dict() 513 for trans in self.trans_list: 514 molec_id = trans.molecule.getModelId() 515 molec = trans.molecule 516 if not molec_id: 517 self.trans_bools.append(False) 518 trans.setModelId('') 519 elif not self.sph_db[self.model_id].has_key(molec_id): 520 trans.setModelId(molec_id) 521 nd = dict([(str(trans),trans.makeDict(1))]) 522 self.sph_db[self.model_id][molec_id] = dict([(molec_id,nd)]) 523 self.sph_db.addChangedKey(self.model_id) 524 self.trans_bools.append(False) 525 else: 526 tdb = self.sph_db[self.model_id][molec_id] 527 gd_ids = [k for k,v in sorted(tdb.items()) 528 if str(trans) in v.keys()] 529 for trans_id in gd_ids: 530 db_trans_dict = self.sph_db[self.model_id][molec_id]\ 531 [trans_id][str(trans)].copy() 532 if self.cCL(this_list=trans.makeDict(),\ 533 modellist=db_trans_dict,code='sphinx'): 534 trans.setModelId(trans_id) 535 self.trans_bools.append(True) 536 if not self.vic is None \ 537 and db_trans_dict.has_key('IN_PROGRESS'): 538 self.vic.addTransInProgress(trans) 539 print 'Sphinx model is currently being '+\ 540 'calculated for %s of %s with ID %s.'\ 541 %(str(trans),molec.molecule,\ 542 trans.getModelId()) 543 elif self.vic is None \ 544 and db_trans_dict.has_key('IN_PROGRESS'): 545 self.addTransInProgress(trans) 546 print 'Sphinx model is currently being ' + \ 547 'calculated in a different CC modeling '+\ 548 'session for %s of %s with ID %s.'\ 549 %(str(trans),molec.molecule,\ 550 trans.getModelId()) 551 else: 552 print 'Sphinx model has been calculated before '+ \ 553 'for %s of %s with ID %s.'\ 554 %(str(trans),molec.molecule,\ 555 trans.getModelId()) 556 break 557 558 #-- No match found. Calculate trans anew. Now give it an id 559 # matching existing ids in terms of pars, or make a new one 560 if trans.getModelId() is None: 561 self.trans_bools.append(False) 562 molec_dbd = self.sph_db[self.model_id][molec_id] 563 564 for i,(k,v) in enumerate(sorted(molec_dbd.items())): 565 #-- No trans in this id yet. 566 # This should NEVER occur! Raise error. 567 # Such empty dicts are removed upon fail/success check 568 # in doSphinx/checkSphinxOutput. 569 if not v.keys(): 570 raise KeyError('Empty trans id found in the ' + \ 571 'database. This should not be ' + \ 572 'possible.') 573 574 #-- Trans already in this id. Keep searching, unless it 575 # is the last entry. Then we have to make a new id 576 if str(trans) in v.keys() \ 577 and i != len(molec_dbd.keys())-1: 578 continue 579 580 #-- Now check if par dict is same except the transition 581 # definition. Trick cCL by adding the db trans def to 582 # this trans' dict. If match, choose this ID 583 # Make sure trans is not in v.keys() already (to avoid 584 # the case of this being the last entry of the db) 585 if not str(trans) in v.keys(): 586 tt = v.keys()[0] 587 extra_dict={'TRANSITION':v[tt]['TRANSITION']} 588 if self.cCL(this_list=trans.makeDict(),\ 589 modellist=v[tt],code='sphinx',\ 590 extra_dict=extra_dict): 591 break 592 593 #-- If end of dict is reached, no proper match was found 594 # Create new id and put it in k. 595 if i == len(molec_dbd.keys())-1: 596 #-- First check if the original molec id is already 597 # in use (it may have been deleted if first use 598 # failed). Never copy output for this though, 599 # since it is the original 600 if not molec_dbd.has_key(molec_id): 601 k = molec_id 602 copied_molecs.append((molec.molecule,k)) 603 else: 604 k = self.makeNewId() 605 self.makeIdLog(new_id=k,molec_id=molec_id) 606 self.sph_db[self.model_id][molec_id][k] = dict() 607 608 #-- It is possible molec_dbd is actually empty. Then use the 609 # molec id as trans id. No need for an id log in this case 610 # nor any copying of output. 611 if not molec_dbd.keys(): 612 k = molec_id 613 copied_molecs.append((molec.molecule,k)) 614 self.sph_db[self.model_id][molec_id][molec_id] = dict() 615 616 #-- The last k value (or the new one in case end of dict was 617 # reached) is set as the model id. 618 trans.setModelId(k) 619 620 #-- When the trans is not yet calculated, either a new id 621 # was made or an existing one is used. It still needs to 622 # be checked if all mline and cooling info is available. 623 # You only want to do this once per session for each 624 # molecule, because ls/ln checks add a lot of overhead. 625 # copyOutput double checks if links already exist 626 if (molec.molecule,k) not in copied_molecs: 627 self.copyOutput(trans,molec_id,k) 628 copied_molecs.append((molec.molecule,k)) 629 td = trans.makeDict(1) 630 self.sph_db[self.model_id][molec_id][k][str(trans)] = td 631 self.sph_db.addChangedKey(self.model_id) 632 633 sph_dbfile.close() 634 if not self.single_session: self.sph_db.sync()

635 636 637

638 - def copyOutput(self,entry,old_id,new_id):

639 640 ''' 641 Copy modelling output based on model_id. 642 643 @param entry: the modeling object for which output is copied 644 @type entry: Molecule() or Transition() 645 @param old_id: The old model_id 646 @type old_id: string 647 @param new_id: the new_model_id 648 @type new_id: string 649 650 ''' 651 652 folder_old = os.path.join(cc.path.gout,'models',old_id) 653 folder_new = os.path.join(cc.path.gout,'models',new_id) 654 lsprocess = subprocess.Popen('ls %s'%folder_old,shell=True,\ 655 stdout=subprocess.PIPE) 656 lsfile = lsprocess.communicate()[0].split('\n') 657 lsfile = [os.path.split(line)[1] 658 for line in lsfile 659 if ((line[0:2] == 'ml' or line[0:4] == 'cool') \ 660 and not entry.isMolecule()) \ 661 or line[0:7] == 'coolfgr' \ 662 or line[0:4] == 'para' \ 663 or line[0:5] == 'input'] 664 if not entry.isMolecule(): 665 lsfile = [line 666 for line in lsfile 667 if not ((line[0:2] == 'ml' \ 668 and os.path.splitext(line)[0].split('_')[-1]\ 669 != entry.molecule.molecule))] 670 671 lsfile = [line 672 for line in lsfile 673 if not (line[0:4] == 'cool' \ 674 and (line.split('_')[-1].replace('.dat','') \ 675 != entry.molecule.molecule \ 676 or line.split('_')[-1].replace('.dat','')=='sampling'\ 677 or line[0:7] == 'coolfgr'))] 678 679 new_lsfile = [line.replace(old_id,new_id) for line in lsfile] 680 DataIO.testFolderExistence(folder_new) 681 lsprocess = subprocess.Popen('ls %s'%folder_new,shell=True,\ 682 stdout=subprocess.PIPE) 683 already_done = lsprocess.communicate()[0].split('\n') 684 for ls,nls in zip(lsfile,new_lsfile): 685 if not nls in already_done: 686 subprocess.call(['ln -s %s %s'%(os.path.join(folder_old,ls),\ 687 os.path.join(folder_new,nls))],\ 688 shell=True)

689 690 691

692 - def doCooling(self,star):

693 694 """ 695 Run Cooling. 696 697 First, database is checked for retrieval of old model. 698 699 @param star: The parameter set for this session 700 @type star: Star() 701 702 """ 703 704 #-- Collect H2O and CO molecule definitions for inclusion in the 705 # cooling inputfile. Also includes abundance_filename info for H2O if 706 # requested 707 if not star.getMolecule('1H1H16O') is None: 708 h2o_dict = star.getMolecule('1H1H16O').makeDict() 709 else: 710 h2o_dict = Molecule('1H1H16O',45,45,648,50).makeDict() 711 if not star.getMolecule('12C16O') is None: 712 co_dict = star.getMolecule('12C16O').makeDict() 713 else: 714 co_dict = Molecule('12C16O',61,61,240,50).makeDict() 715 716 #-- no abundance profiles should be possible for CO. 717 if co_dict.has_key('MOLECULE_TABLE'): 718 raise IOError('CO cannot be attributed a custom abundance ' + \ 719 'profile at this time.') 720 721 #-- F_H2O is irrelevant if an abundance file is passed for oH2O 722 if h2o_dict.has_key('MOLECULE_TABLE'): 723 del self.command_list['F_H2O'] 724 725 #-- Collect all H2O molecular information important for cooling 726 molec_dict = dict([(k,h2o_dict[k]) 727 for k in self.cooling_molec_keys 728 if h2o_dict.has_key(k)]) 729 730 #-- Check database: only include H2O extra keywords if 731 # abundance_filename is present. CO can't have this anyway. 732 model_bool = self.checkCoolingDatabase(molec_dict=molec_dict.copy()) 733 734 #- Run cooling if above is False 735 if not model_bool: 736 DataIO.testFolderExistence(os.path.join(cc.path.gout,'models',\ 737 self.model_id)) 738 commandfile = ['%s=%s'%(k,v) 739 for k,v in sorted(self.command_list.items()) 740 if k != 'R_POINTS_MASS_LOSS'] + \ 741 ['####'] + \ 742 ['%s=%s'%('MOLECULE',co_dict['MOLECULE'])] + \ 743 ['%s=%s'%(k,h2o_dict[k]) 744 for k in self.cooling_molec_keys + ['MOLECULE'] 745 if h2o_dict.has_key(k)] + ['####'] 746 if self.command_list.has_key('R_POINTS_MASS_LOSS'): 747 commandfile.extend(['%s=%s'%('R_POINTS_MASS_LOSS',v) 748 for v in self.command_list\ 749 ['R_POINTS_MASS_LOSS']] + \ 750 ['####']) 751 filename = os.path.join(cc.path.gout,'models',\ 752 'gastronoom_' + self.model_id + '.inp') 753 DataIO.writeFile(filename,commandfile) 754 if not self.skip_cooling: 755 self.execGastronoom(subcode='cooling',filename=filename) 756 self.cool_done = True 757 if os.path.isfile(os.path.join(cc.path.gout,'models',\ 758 self.model_id,'coolfgr_all%s.dat'\ 759 %self.model_id)): 760 #-- Note that there is no need to create a log parameter file 761 # since the inputfiles are still available, and they always 762 # contain all cooling keywords. (maybe change for h2o cooling) 763 if self.cool_db[self.model_id].has_key('IN_PROGRESS'): 764 del self.cool_db[self.model_id]['IN_PROGRESS'] 765 self.cool_db.addChangedKey(self.model_id) 766 else: 767 print 'Cooling model calculation failed. No entry is added '+ \ 768 'to the database.' 769 del self.cool_db[self.model_id] 770 self.model_id = '' 771 if not self.single_session: self.cool_db.sync()

772 773 774

775 - def doMline(self,star):

776 777 """ 778 Run mline. 779 780 First, database is checked for retrieval of old models. 781 782 @param star: The parameter set for this session 783 @type star: Star() 784 785 """ 786 787 #-- Make sure to reset this in case an iteration between cooling and 788 # mline is happening 789 self.mline_done = False 790 model_bools = self.checkMlineDatabase() 791 del self.command_list['R_OUTER'] 792 del self.command_list['OUTER_R_MODE'] 793 for molec,model_bool in zip(self.molec_list,model_bools): 794 if not model_bool: 795 self.updateModel(molec.getModelId()) 796 commandfile = ['%s=%s'%(k,v) 797 for k,v in sorted(self.command_list.items()) 798 if k != 'R_POINTS_MASS_LOSS'] +\ 799 ['####'] + \ 800 ['%s=%s'%(k,v) 801 for k,v in sorted(molec.makeDict().items())] +\ 802 ['####'] 803 if self.command_list.has_key('R_POINTS_MASS_LOSS'): 804 commandfile.extend(['%s=%s'%('R_POINTS_MASS_LOSS',v) 805 for v in self.command_list\ 806 ['R_POINTS_MASS_LOSS']] +\ 807 ['####']) 808 filename = os.path.join(cc.path.gout,'models',\ 809 'gastronoom_%s.inp'%molec.getModelId()) 810 DataIO.writeFile(filename,commandfile) 811 self.execGastronoom(subcode='mline',filename=filename) 812 self.mline_done=True 813 path = os.path.join(cc.path.gout,'models',molec.getModelId()) 814 fns = 'ml*{}_{}.dat'.format(molec.getModelId(),molec.molecule) 815 if len(glob(os.path.join(path,fns))) == 3: 816 #-- Remove in-progress entry. 817 if self.ml_db[self.model_id][molec.getModelId()]\ 818 [molec.molecule].has_key('IN_PROGRESS'): 819 del self.ml_db[self.model_id][molec.getModelId()]\ 820 [molec.molecule]['IN_PROGRESS'] 821 822 #-- Write mline keywords not included in sph files but used 823 # in the database in an extra log file. Only do this if it 824 # doesn't already exist. The parameters should be the same 825 # for all transitions with this model id. 826 mlfn = 'mline_parameters_{}.log'.format(molec.molecule) 827 mlfn = os.path.join(path,mlfn) 828 if not os.path.isfile(mlfn): 829 #-- Add MOLECULE too. Cuz, why not. For TRANSITION, that 830 # info is recreated from sph files. Not so for mline. 831 mlfile = ['{}={}'.format(k,v) 832 for k,v in sorted(molec.makeDict().items())] 833 DataIO.writeFile(mlfn,mlfile) 834 else: 835 del self.ml_db[self.model_id][molec.getModelId()]\ 836 [molec.molecule] 837 #-- Remove the molecule id if it does not contain molecules 838 # anymore. The id is thus unused. 839 if not self.ml_db[self.model_id][molec.getModelId()].keys(): 840 del self.ml_db[self.model_id][molec.getModelId()] 841 print 'Mline model calculation failed for'\ 842 '%s. No entry is added to the database.'\ 843 %(molec.molecule) 844 molec.setModelId('') 845 846 #-- Synchronize db: Both when successful or failure. 847 self.ml_db.addChangedKey(self.model_id) 848 if not self.single_session: self.ml_db.sync() 849 850 851 if set([molec.getModelId() for molec in self.molec_list]) == set(['']): 852 #- no mline models calculated: stop GASTRoNOoM here 853 self.model_id = '' 854 print 'Mline model calculation failed for all requested ' + \ 855 'molecules. Stopping GASTRoNOoM here!' 856 else: 857 #- at least one molecule was successfully calculated, so start 858 #- Sphinx, hence if vic is requested, the cooling model_id can now 859 #- be added to the models list 860 if not self.vic is None and self.sphinx: 861 #- add the command list to the vic models list 862 self.vic.addModel(self.model_id,self.command_list)

863 864 865

866 - def doSphinx(self,star):

867 868 """ 869 Run Sphinx. 870 871 First, database is checked for retrieval of old models. 872 873 @param star: The parameter set for this session 874 @type star: Star() 875 876 """ 877 878 self.checkSphinxDatabase() 879 print '%i transitions out of %i not yet calculated.'\ 880 %(len([boolean for boolean in self.trans_bools if not boolean]),\ 881 len(self.trans_bools)) 882 for i,(trans_bool,trans) in enumerate(zip(self.trans_bools,\ 883 self.trans_list)): 884 if not trans_bool and trans.getModelId(): 885 if not self.sphinx: 886 #- Only transitions with no db entry will get empty model id 887 del self.sph_db[self.model_id][trans.molecule.getModelId()]\ 888 [trans.getModelId()][str(trans)] 889 self.sph_db.addChangedKey(self.model_id) 890 trans.setModelId('') 891 elif not self.vic is None: 892 #- add transition to the vic translist for this cooling id 893 self.vic.addTrans(trans) 894 elif self.recover_sphinxfiles: 895 self.checkSphinxOutput(trans) 896 else: 897 self.updateModel(trans.getModelId()) 898 commandfile = ['%s=%s'%(k,v) 899 for k,v in sorted(self.command_list.items()) 900 if k != 'R_POINTS_MASS_LOSS'] + ['####'] + \ 901 ['%s=%s'%(k,v) 902 for k,v in sorted(trans.molecule.makeDict()\ 903 .items())] + \ 904 ['####'] + \ 905 ['%s=%s'%(k,v) 906 for k,v in sorted(trans.makeDict()\ 907 .items())] + \ 908 ['######'] 909 if self.command_list.has_key('R_POINTS_MASS_LOSS'): 910 commandfile.extend(['%s=%s'%('R_POINTS_MASS_LOSS',v) 911 for v in self.command_list\ 912 ['R_POINTS_MASS_LOSS']] + \ 913 ['####']) 914 filename = os.path.join(cc.path.gout,'models',\ 915 'gastronoom_%s.inp'\ 916 %trans.getModelId()) 917 DataIO.writeFile(filename,commandfile) 918 print 'Starting calculation for transition %i out of %i.'\ 919 %(i+1,len(self.trans_bools)) 920 self.execGastronoom(subcode='sphinx',filename=filename) 921 self.checkSphinxOutput(trans) 922 if not self.single_session: self.sph_db.sync() 923 924 #- check if at least one of the transitions was calculated: then 925 #- self.model_id doesnt have to be changed 926 self.finalizeSphinx() 927 928 #-- Sync the sphinx db in case self.sphinx is False or 929 # self.recover_sphinxfiles is True, to make sure sph_db is up-to-date. 930 # In case models are calculated or vic is ran this is done in other 931 # places in the code. 932 if not self.sphinx or self.recover_sphinxfiles: 933 self.sph_db.sync() 934 935 mline_not_available = set([trans.molecule.getModelId() 936 for boolean,trans in zip(self.trans_bools,\ 937 self.trans_list) 938 if not boolean]) \ 939 == set(['']) 940 if not self.vic is None and self.sphinx and (False in self.trans_bools \ 941 and not mline_not_available): 942 self.vic.queueModel() 943 elif not self.vic is None and self.sphinx and (False not in self.trans_bools\ 944 or mline_not_available): 945 self.vic.reset()

946 947 948

949 - def finalizeSphinx(self):

950 951 ''' 952 Check if at least one of the transitions has been calculated correctly. 953 954 if not, self.model_id is set to "". 955 956 ''' 957 958 if set([trans.getModelId() for trans in self.trans_list]) == set(['']): 959 self.model_id = '' 960 for trans in self.trans_in_progress: 961 if not self.sph_db[self.model_id][trans.molecule.getModelId()]\ 962 [trans.getModelId()].has_key(str(trans)) \ 963 or self.sph_db[self.model_id][trans.molecule.getModelId()]\ 964 [trans.getModelId()][str(trans)]\ 965 .has_key('IN_PROGRESS'): 966 trans.setModelId('')

967 968 969

970 - def checkSphinxOutput(self,trans):

971 972 ''' 973 Check if sphinx output is complete and update the database with 974 calculated models. 975 976 Requires model_id and path defined in Gastronoom instance. 977 978 @param trans: the transition that is being checked 979 @type trans: Transition() 980 981 ''' 982 983 filename = trans.makeSphinxFilename(number='*') 984 path = os.path.join(cc.path.gout,'models',trans.getModelId()) 985 #- Sphinx puts out 2 files per transition 986 if len(glob(os.path.join(path,filename))) == 2: 987 if self.sph_db[self.model_id][trans.molecule.getModelId()]\ 988 [trans.getModelId()][str(trans)]\ 989 .has_key('IN_PROGRESS'): 990 del self.sph_db[self.model_id][trans.molecule.getModelId()]\ 991 [trans.getModelId()][str(trans)]['IN_PROGRESS'] 992 993 #-- Write sphinx keywords not included in sph filenames but used in 994 # the database in an extra log file. Only do this if it doesn't 995 # already exist. The parameters should be the same for all 996 # transitions with this model id. 997 sphfn = os.path.join(path,'sphinx_parameters.log') 998 if not os.path.isfile(sphfn): 999 sphfile = ['{}={}'.format(k,v) 1000 for k,v in sorted(trans.makeDict().items()) 1001 if k != 'TRANSITION'] 1002 DataIO.writeFile(sphfn,sphfile) 1003 1004 #-- Print to the shell that the line is finished. 1005 print 'Sphinx model calculated successfully for '+\ 1006 '%s of %s with id %s.'%(str(trans),trans.molecule.molecule,\ 1007 trans.getModelId()) 1008 else: 1009 del self.sph_db[self.model_id][trans.molecule.getModelId()]\ 1010 [trans.getModelId()][str(trans)] 1011 #-- Remove the transition id if it does not contain transitions 1012 # anymore. The id is thus unused. 1013 if not self.sph_db[self.model_id][trans.molecule.getModelId()]\ 1014 [trans.getModelId()].keys(): 1015 del self.sph_db[self.model_id][trans.molecule.getModelId()]\ 1016 [trans.getModelId()] 1017 print 'Sphinx model calculation failed for %s of %s with id %s.'\ 1018 %(str(trans),trans.molecule.molecule,trans.getModelId()) 1019 print 'No entry is added to the Sphinx database.' 1020 trans.setModelId('') 1021 1022 self.sph_db.addChangedKey(self.model_id)

1023 1024

1025 - def setCommandKey(self,comm_key,star,star_key=None,alternative=None):

1026 1027 ''' 1028 Try setting a key in the command_list from a star instance. 1029 1030 If the key is unknown, it is left open and will be filled in from the 1031 standard gastronoom inputfile. 1032 1033 @param comm_key: the name of the keyword in the command list 1034 @type comm_key: string 1035 @param star: Parameter set for this session 1036 @type star: Star() 1037 1038 @keyword star_key: the name of the keyword in Star() (minus '_%s' 1039 %key_type (DUST or GAS), which is added as well in a 1040 second attempt if the first without the addition is 1041 not found). If None, it's equal to comm_key 1042 1043 (default: None) 1044 @type star_key: string 1045 @keyword alternative: a default value passed from the standard 1046 inputfile that is used if the keyword or the 1047 keyword + '_%s'%key_type is not found in Star(). 1048 1049 (default: None) 1050 @type alternative: string 1051 1052 @return: Success? 1053 @rtype: bool 1054 1055 ''' 1056 1057 keyword_int_list = ['ITERA_COOLING','LOG_DEPTH_STEP_POWER',\ 1058 'USE_MLINE_COOLING_RATE_CO',\ 1059 'USE_NEW_DUST_KAPPA_FILES','STEP_RIN_ROUT',\ 1060 'STEP_RS_RIN'] 1061 exp_not_list = ['STEP_RIN_ROUT','STEP_RS_RIN'] 1062 make_int = comm_key in keyword_int_list 1063 1064 #- Make sure large integers are given in exponential notation 1065 #- Fortran cannot work with very large integers. 1066 exp_not = comm_key in exp_not_list 1067 1068 return super(Gastronoom, self).setCommandKey(comm_key,star,'GAS',\ 1069 star_key,alternative,\ 1070 make_int,exp_not)

1071 1072 1073

1074 - def doGastronoom(self,star):

1075 1076 """ 1077 Run GASTRoNOoM-cooling. 1078 1079 The input parameter list is prepared here. 1080 1081 @param star: Parameter set for this session 1082 @type star: Star() 1083 1084 """ 1085 1086 print '***********************************' 1087 print '** Making input file for GASTRoNOoM' 1088 1089 #-- Add the previous cooling model_id to the list of new entries, so it 1090 # does not get deleted if replace_db_entry == 1. 1091 # This id is not known in new_entries, as the new_entries are passed 1092 # for the previous models, not the current one.(ie when iterations>1) 1093 if self.model_id: 1094 self.new_entries.append(self.model_id) 1095 1096 #-- Make sure to reset this in case an iteration between cooling and 1097 # mline, cooling and mcmax is happening 1098 self.cool_done = False 1099 self.model_id = self.makeNewId() 1100 self.trans_list=star['GAS_LINES'] 1101 self.molec_list=star['GAS_LIST'] 1102 self.command_list = dict() 1103 self.command_list['DATA_DIRECTORY'] = '"' + cc.path.gdata + '"' 1104 self.command_list['OUTPUT_DIRECTORY'] \ 1105 = '"' + os.path.join(cc.path.gout,'models',self.model_id) +'/"' 1106 self.command_list['PARAMETER_FILE'] \ 1107 = '"' + os.path.join(cc.path.gout,'models',self.model_id,\ 1108 'parameter_file_%s.dat'%self.model_id)+'"' 1109 self.command_list['OUTPUT_SUFFIX'] = self.model_id 1110 1111 if star['TEMDUST_FILENAME'] == 'temdust.kappa': 1112 self.command_list['SPEC_DENS_DUST'] = 3.3 1113 else: 1114 self.command_list['SPEC_DENS_DUST'] = star['SPEC_DENS_DUST'] 1115 1116 self.command_list['DUST_TEMPERATURE_FILENAME'] \ 1117 = '"%s"'%star['DUST_TEMPERATURE_FILENAME'] 1118 self.command_list['TEMDUST_FILENAME'] = '"%s"'%star['TEMDUST_FILENAME'] 1119 self.command_list['R_STAR'] = float(star['R_STAR'])*star.Rsun 1120 1121 if star['MDOT_MODE'].upper() != 'CONSTANT': 1122 self.command_list['USE_DENSITY_NON_CONSISTENT'] \ 1123 = str(int(star['USE_DENSITY_NON_CONSISTENT'])) 1124 if star['MDOT_MODE'].upper() == 'OTHER': 1125 try: 1126 self.command_list['R_POINTS_MASS_LOSS'] \ 1127 = star['R_POINTS_MASS_LOSS'] 1128 except KeyError: 1129 raise KeyError('R_POINTS_MASS_LOSS not present. Include '+\ 1130 'this parameter in your inputfile or set '+\ 1131 'MDOT_MODE to something different than OTHER.') 1132 1133 #- always has to be included, 0.5 if not present in the inputfile, 1134 #- also sets power law innermost region 1135 self.command_list['TEMPERATURE_EPSILON'] \ 1136 = star['TEMPERATURE_EPSILON_GAS'] 1137 1138 #- the next few keywords only if the mode is not cooling, but epsilon 1139 if star['TEMPERATURE_MODE_GAS'] != 'cooling': 1140 if float(star['TEMPERATURE_EPSILON2_GAS']): 1141 self.command_list['TEMPERATURE_EPSILON2'] \ 1142 = star['TEMPERATURE_EPSILON2_GAS'] 1143 self.command_list['RADIUS_EPSILON2'] \ 1144 = star['RADIUS_EPSILON2_GAS'] 1145 #- Only add third power law if the second power law is non-zero 1146 if float(star['TEMPERATURE_EPSILON3_GAS']): 1147 self.command_list['TEMPERATURE_EPSILON3'] \ 1148 = star['TEMPERATURE_EPSILON3_GAS'] 1149 self.command_list['RADIUS_EPSILON3'] \ 1150 = star['RADIUS_EPSILON3_GAS'] 1151 1152 #-- Add the stochastic velocity filename if mode is filename, or power 1153 # and inner value if mode is powerlaw. Mode and Vel are always added. 1154 if star['STOCHASTIC_VEL_MODE'].upper() == 'POWERLAW': 1155 self.command_list['STOCHASTIC_VEL_POWER'] \ 1156 = star['STOCHASTIC_VEL_POWER'] 1157 self.command_list['STOCHASTIC_VEL_INNER'] \ 1158 = star['STOCHASTIC_VEL_INNER'] 1159 elif star['STOCHASTIC_VEL_MODE'].upper() == 'FILENAME': 1160 self.command_list['STOCHASTIC_VEL_FILENAME'] \ 1161 = star['STOCHASTIC_VEL_FILENAME'] 1162 1163 self.setCommandKey('DUST_TO_GAS',star,'DUST_TO_GAS_INITIAL',\ 1164 self.standard_inputfile['DUST_TO_GAS']) 1165 1166 #- Both pars set explicitly because they are also present in 1167 #- mline_keywords 1168 self.setCommandKey('R_OUTER',star,'R_OUTER',\ 1169 self.standard_inputfile['R_OUTER']) 1170 self.setCommandKey('OUTER_R_MODE',star,'OUTER_R_MODE',\ 1171 self.standard_inputfile['OUTER_R_MODE']) 1172 add_keys = [k 1173 for k in self.standard_inputfile.keys() 1174 if not (self.command_list.has_key(k) \ 1175 or k in self.mline_keywords + self.sphinx_keywords)] 1176 [self.setCommandKey(k,star,alternative=self.standard_inputfile[k]) 1177 for k in add_keys] 1178 print '** DONE!' 1179 print '***********************************' 1180 1181 #- model/output naming entries are excluded when comparing new models 1182 #- with database models 1183 #- Start the model calculation 1184 self.doCooling(star)

Source Code for Module ComboCode.cc.modeling.codes.Gastronoom