solvcon.kerpak.cese

1 # -*- coding: UTF-8 -*- 2 # 3 # Copyright (C) 2008-2011 Yung-Yu Chen <yyc@solvcon.net>. 4 # 5 # This program is free software; you can redistribute it and/or modify 6 # it under the terms of the GNU General Public License as published by 7 # the Free Software Foundation; either version 2 of the License, or 8 # (at your option) any later version. 9 # 10 # This program is distributed in the hope that it will be useful, 11 # but WITHOUT ANY WARRANTY; without even the implied warranty of 12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 # GNU General Public License for more details. 14 # 15 # You should have received a copy of the GNU General Public License along 16 # with this program; if not, write to the Free Software Foundation, Inc., 17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 18 19 """ 20 Space-time Conservation Element and Solution Element (CESE) method. 21 """ 22 23 from ctypes import Structure 24 from solvcon.solver import BlockSolver 25 from solvcon.case import BlockCase 26 from solvcon.boundcond import BC, periodic 27 from solvcon.anchor import Anchor 28 from solvcon.hook import Hook, BlockHook

29 30 ################################################################################ 31 # Solver 32 ################################################################################ 33 34 -class CeseSolverExedata(Structure):

35 """ 36 Data structure to interface with C. 37 """ 38 from ctypes import c_int, c_double, c_void_p 39 _fields_ = [ 40 # inherited. 41 ('ncore', c_int), ('neq', c_int), 42 ('time', c_double), ('time_increment', c_double), 43 # mesh shape. 44 ('ndim', c_int), ('nnode', c_int), ('nface', c_int), ('ncell', c_int), 45 ('nbound', c_int), 46 ('ngstnode', c_int), ('ngstface', c_int), ('ngstcell', c_int), 47 # group shape. 48 ('ngroup', c_int), ('gdlen', c_int), 49 # parameter shape. 50 ('nsca', c_int), ('nvec', c_int), 51 # function pointer. 52 ('jacofunc', c_void_p), ('taufunc', c_void_p), ('omegafunc', c_void_p), 53 # scheme. 54 ('alpha', c_int), ('taylor', c_double), 55 ('cnbfac', c_double), ('sftfac', c_double), 56 ('taumin', c_double), ('taumax', c_double), ('tauscale', c_double), 57 ('omegamin', c_double), ('omegascale', c_double), 58 ('mqmin', c_double), ('mqscale', c_double), 59 # meta array. 60 ('fctpn', c_void_p), 61 ('cltpn', c_void_p), ('clgrp', c_void_p), 62 ('grpda', c_void_p), 63 # metric array. 64 ('ndcrd', c_void_p), 65 ('fccnd', c_void_p), ('fcnml', c_void_p), 66 ('clcnd', c_void_p), ('clvol', c_void_p), 67 ('cecnd', c_void_p), ('cevol', c_void_p), 68 ('mqlty', c_void_p), 69 # connectivity array. 70 ('fcnds', c_void_p), ('fccls', c_void_p), 71 ('clnds', c_void_p), ('clfcs', c_void_p), 72 # solution array. 73 ('sol', c_void_p), ('dsol', c_void_p), 74 ('solt', c_void_p), 75 ('soln', c_void_p), ('dsoln', c_void_p), 76 ('cfl', c_void_p), ('ocfl', c_void_p), 77 ('amsca', c_void_p), ('amvec', c_void_p), 78 ] 79 del c_int, c_double, c_void_p 80

81 - def __set_pointer(self, svr, aname, shf):

82 from ctypes import c_void_p 83 ptr = getattr(svr, aname)[shf:].ctypes._as_parameter_ 84 setattr(self, aname, ptr)

85

86 - def __init__(self, *args, **kw):

87 from ctypes import c_int, c_double, POINTER, c_void_p, byref, cast 88 svr = kw.pop('svr', None) 89 super(CeseSolverExedata, self).__init__(*args, **kw) 90 if svr == None: 91 return 92 # inherited. 93 for key in ('ncore', 'neq', 'time', 'time_increment',): 94 setattr(self, key, getattr(svr, key)) 95 # mesh shape. 96 for key in ('ndim', 'nnode', 'nface', 'ncell', 'nbound', 97 'ngstnode', 'ngstface', 'ngstcell',): 98 setattr(self, key, getattr(svr, key)) 99 # group shape. 100 for key in ('ngroup', 'gdlen',): 101 setattr(self, key, getattr(svr, key)) 102 # parameter shape. 103 for key in ('nsca', 'nvec',): 104 setattr(self, key, getattr(svr, key)) 105 # function pointer. 106 self.jacofunc = cast(svr._jacofunc_, c_void_p).value 107 self.taufunc = cast(getattr(svr._clib_cese, 108 'calc_tau_'+svr.tauname), c_void_p).value 109 self.omegafunc = cast(getattr(svr._clib_cese, 110 'calc_omega_'+svr.omeganame), c_void_p).value 111 # scheme. 112 for key in ('alpha', 'taylor', 'cnbfac', 'sftfac', 113 'taumin', 'taumax', 'tauscale', 'omegamin', 'omegascale', 114 'mqmin', 'mqscale',): 115 setattr(self, key, getattr(svr, key)) 116 # meta array. 117 for aname in ('fctpn',): 118 self.__set_pointer(svr, aname, svr.ngstface) 119 for aname in ('cltpn', 'clgrp',): 120 self.__set_pointer(svr, aname, svr.ngstcell) 121 for aname in ('grpda',): 122 self.__set_pointer(svr, aname, 0) 123 # metric array. 124 for aname in ('ndcrd',): 125 self.__set_pointer(svr, aname, svr.ngstnode) 126 for aname in ('fccnd', 'fcnml',): 127 self.__set_pointer(svr, aname, svr.ngstface) 128 for aname in ('clcnd', 'clvol', 'cecnd', 'cevol',): 129 self.__set_pointer(svr, aname, svr.ngstcell) 130 # mesh quality. 131 for aname in ('mqlty',): 132 self.__set_pointer(svr, aname, svr.ngstcell) 133 # connectivity array. 134 for aname in ('fcnds', 'fccls',): 135 self.__set_pointer(svr, aname, svr.ngstface) 136 for aname in ('clnds', 'clfcs',): 137 self.__set_pointer(svr, aname, svr.ngstcell) 138 # solution array. 139 for aname in ('sol', 'dsol', 'solt', 'soln', 'dsoln',): 140 self.__set_pointer(svr, aname, svr.ngstcell) 141 for aname in ('cfl', 'ocfl'): 142 self.__set_pointer(svr, aname, 0) 143 for aname in ('amsca', 'amvec'): 144 self.__set_pointer(svr, aname, svr.ngstcell)

145

146 -class CeseSolver(BlockSolver):

147 """ 148 The base solver class for multi-dimensional CESE code. 149 150 @cvar _gdlen_: length per group data. Must be overridden. 151 @ctype _gdlen_: int 152 @cvar _jacofunc_: ctypes function to Jacobian calculator. Must be 153 overridden. 154 @ctype _jacofunc_: ctypes.FuncPtr 155 156 @ivar debug: flag for debugging. 157 @itype debug: bool 158 159 @ivar diffname: name of gradient calculation function; tau is default, 160 omega is selectable. 161 @itype diffname: str 162 @ivar tauname: name of tau function; default linear. 163 @itype tauname: str 164 @ivar omeganame: name of omega function; default scale. 165 @itype omeganame: str 166 @ivar alpha: parameter to the weighting function. 167 @itype alpha: int 168 @ivar taylor: factor for Taylor's expansion; 0 off, 1 on. 169 @itype taylor: float 170 @ivar cnbfac: factor to use BCE centroid, othersize midpoint; 0 off, 1 on. 171 @itype cnbfac: float 172 @ivar sftfac: factor to shift gradient shape; 0 off, 1 on. 173 @itype sftfac: float 174 @ivar taumin: the lower bound of tau. 175 @itype taumin: float 176 @ivar taumax: the upper bound of tau. 177 @itype taumax: float 178 @ivar tauscale: scaling of tau. 179 @itype tauscale: float 180 @ivar omegamin: the lower bound of omega. 181 @itype omegamin: float 182 @ivar omegascale: scaling of omega. 183 @itype omegascale: float 184 185 @ivar grpda: group data. 186 @ivar cecnd: solution points for CCEs and BCEs. 187 @ivar cevol: CCE and BCE volumes. 188 @ivar solt: temporal diffrentiation of solution. 189 @ivar sol: current solution. 190 @ivar soln: next solution. 191 @ivar dsol: current gradient of solution. 192 @ivar dsoln: next gradient of solution. 193 @ivar cfl: CFL number. 194 @ivar ocfl: original CFL number. 195 @ivar amsca: Parameter scalar array. 196 @ivar amvec: Parameter vector array. 197 """ 198 199 _exedatatype_ = CeseSolverExedata 200 _interface_init_ = ['cecnd', 'cevol'] 201 _solution_array_ = ['sol', 'soln', 'dsol', 'dsoln'] 202 203 _gdlen_ = None 204 _jacofunc_ = None 205

206 - def __init__(self, blk, *args, **kw):

207 from numpy import empty 208 self.debug = kw.pop('debug', False) 209 nsca = kw.pop('nsca', 0) 210 nvec = kw.pop('nvec', 0) 211 diffname = kw.pop('diffname', None) 212 self.diffname = diffname if diffname != None else 'tau' 213 tauname = kw.pop('tauname', None) 214 self.tauname = tauname if tauname != None else 'linear' 215 omeganame = kw.pop('omeganame', None) 216 self.omeganame = omeganame if omeganame != None else 'scale' 217 self.alpha = int(kw.pop('alpha', 0)) 218 self.taylor = float(kw.pop('taylor', 1.0)) 219 self.cnbfac = float(kw.pop('cnbfac', 1.0)) 220 self.sftfac = float(kw.pop('sftfac', 1.0)) 221 self.taumin = float(kw.pop('taumin', 0.0)) 222 self.taumax = float(kw.pop('taumax', 1.0)) 223 self.tauscale = float(kw.pop('tauscale', 0.0)) 224 self.omegamin = float(kw.pop('omegamin', 1.1)) 225 self.omegascale = float(kw.pop('omegascale', 0.0)) 226 self.mqmin = float(kw.pop('mqmin', 1.0)) 227 self.mqscale = float(kw.pop('mqscale', 0.0)) 228 super(CeseSolver, self).__init__(blk, *args, **kw) 229 fpdtype = self.fpdtype 230 ndim = self.ndim 231 ncell = self.ncell 232 ngstcell = self.ngstcell 233 neq = self.neq 234 ngroup = self.ngroup 235 # meta array. 236 self.grpda = empty((ngroup, self._gdlen_), dtype=fpdtype) 237 # dual mesh. 238 self.cecnd = empty((ngstcell+ncell, self.CLMFC+1, ndim), dtype=fpdtype) 239 self.cevol = empty((ngstcell+ncell, self.CLMFC+1), dtype=fpdtype) 240 # mesh quality. 241 self.mqlty = empty(ngstcell+ncell, dtype=fpdtype) 242 # solutions. 243 self.solt = empty((ngstcell+ncell, neq), dtype=fpdtype) 244 self.sol = empty((ngstcell+ncell, neq), dtype=fpdtype) 245 self.soln = empty((ngstcell+ncell, neq), dtype=fpdtype) 246 self.dsol = empty((ngstcell+ncell, neq, ndim), dtype=fpdtype) 247 self.dsoln = empty((ngstcell+ncell, neq, ndim), dtype=fpdtype) 248 self.cfl = empty(ncell, dtype=fpdtype) 249 self.ocfl = empty(ncell, dtype=fpdtype) 250 self.amsca = empty((ngstcell+ncell, nsca), dtype=fpdtype) 251 self.amvec = empty((ngstcell+ncell, nvec, ndim), dtype=fpdtype)

252 253 @property

254 - def gdlen(self):

255 """ 256 Length per group data. 257 """ 258 return self._gdlen_

259 @property

260 - def nsca(self):

261 return self.amsca.shape[1]

262 @property

263 - def nvec(self):

264 return self.amvec.shape[1]

265

266 - def locate_point(self, *args):

267 from ctypes import byref, c_int, c_void_p 268 from numpy import array 269 crd = array(args, dtype=self.fpdtype) 270 picl = c_int(0) 271 pifl = c_int(0) 272 pjcl = c_int(0) 273 pjfl = c_int(0) 274 self._clib_cese.locate_point( 275 byref(self.exd), 276 crd.ctypes._as_parameter_, 277 byref(picl), 278 byref(pifl), 279 byref(pjcl), 280 byref(pjfl), 281 ) 282 return picl.value, pifl.value, pjcl.value, pjfl.value

283

284 - def init(self, **kw):

285 self._tcall(self._clib_cese.calc_ce, 0, self.ncell) 286 self._tcall(self._clib_cese.qualify_mesh, -self.ngstcell, self.ncell) 287 super(CeseSolver, self).init(**kw)

288

289 - def boundcond(self):

290 super(CeseSolver, self).boundcond() 291 self.call_non_interface_bc('soln') 292 self.call_non_interface_bc('dsoln')

293 294 ################################################## 295 # library. 296 ################################################## 297 from ctypes import c_int 298 from solvcon.dependency import getcdll 299 __clib_cese = { 300 2: getcdll('cese2d'), 301 3: getcdll('cese3d'), 302 } 303 del getcdll, c_int 304 @property

305 - def _clib_cese(self):

306 return self.__clib_cese[self.ndim]

307 308 ################################################## 309 # marching algorithm. 310 ################################################## 311 MMNAMES = list() 312 MMNAMES.append('update')

313 - def update(self, worker=None):

314 from ctypes import c_void_p 315 if self.debug: self.mesg('update') 316 # exchange solution and gradient. 317 self.sol, self.soln = self.soln, self.sol 318 self.dsol, self.dsoln = self.dsoln, self.dsol 319 # reset pointers in execution data. 320 ngstcell = self.ngstcell 321 self.exd.sol = self.sol[ngstcell:].ctypes._as_parameter_ 322 self.exd.soln = self.soln[ngstcell:].ctypes._as_parameter_ 323 self.exd.dsol = self.dsol[ngstcell:].ctypes._as_parameter_ 324 self.exd.dsoln = self.dsoln[ngstcell:].ctypes._as_parameter_ 325 if self.debug: self.mesg(' done.\n')

326 327 MMNAMES.append('ibcam')

328 - def ibcam(self, worker=None):

329 if self.debug: self.mesg('ibcam') 330 if worker: 331 if self.nsca: self.exchangeibc('amsca', worker=worker) 332 if self.nvec: self.exchangeibc('amvec', worker=worker) 333 if self.debug: self.mesg(' done.\n')

334 335 MMNAMES.append('calcsolt')

336 - def calcsolt(self, worker=None):

337 if self.debug: self.mesg('calcsolt') 338 self._tcall(self._clib_cese.calc_solt, -self.ngstcell, self.ncell, 339 tickerkey='calcsolt') 340 if self.debug: self.mesg(' done.\n')

341 MMNAMES.append('calcsoln')

342 - def calcsoln(self, worker=None):

343 if self.debug: self.mesg('calcsoln') 344 func = self._clib_cese.calc_soln 345 self._tcall(func, 0, self.ncell, tickerkey='calcsoln') 346 if self.debug: self.mesg(' done.\n')

347 348 MMNAMES.append('ibcsoln')

349 - def ibcsoln(self, worker=None):

350 if self.debug: self.mesg('ibcsoln') 351 if worker: self.exchangeibc('soln', worker=worker) 352 if self.debug: self.mesg(' done.\n')

353 MMNAMES.append('bcsoln')

354 - def bcsoln(self, worker=None):

355 if self.debug: self.mesg('bcsoln') 356 self.call_non_interface_bc('soln') 357 if self.debug: self.mesg(' done.\n')

358 359 MMNAMES.append('calccfl')

360 - def calccfl(self, worker=None):

361 """ 362 @return: mincfl, maxcfl, number of tuned CFL, accumulated number of 363 tuned CFL. 364 @rtype: tuple 365 """ 366 raise NotImplementedError

367 368 MMNAMES.append('calcdsoln')

369 - def calcdsoln(self, worker=None):

370 if self.debug: self.mesg('calcdsoln') 371 func = getattr(self._clib_cese, 'calc_dsoln_'+self.diffname) 372 self._tcall(func, 0, self.ncell, tickerkey='calcdsoln') 373 if self.debug: self.mesg(' done.\n')

374 375 MMNAMES.append('ibcdsoln')

376 - def ibcdsoln(self, worker=None):

377 if self.debug: self.mesg('ibcdsoln') 378 if worker: self.exchangeibc('dsoln', worker=worker) 379 if self.debug: self.mesg(' done.\n')

380 MMNAMES.append('bcdsoln')

381 - def bcdsoln(self, worker=None):

382 if self.debug: self.mesg('bcdsoln') 383 self.call_non_interface_bc('dsoln') 384 if self.debug: self.mesg(' done.\n')

385

386 ################################################################################ 387 # Case. 388 ################################################################################ 389 390 -class CeseCase(BlockCase):

391 """ 392 Base CESE caes. 393 """ 394 defdict = { 395 'execution.verified_norm': -1.0, 396 'solver.debug_cese': False, 397 'solver.diffname': None, 398 'solver.tauname': None, 399 'solver.omeganame': None, 400 'solver.alpha': 1, 401 'solver.taylor': 1.0, 402 'solver.cnbfac': 1.0, 403 'solver.sftfac': 1.0, 404 'solver.taumin': None, 405 'solver.taumax': None, 406 'solver.tauscale': None, 407 'solver.omegamin': None, 408 'solver.omegascale': None, 409 'solver.mqmin': None, 410 'solver.mqscale': None, 411 }

412 - def __init__(self, **kw):

413 import os 414 super(CeseCase, self).__init__(**kw) 415 pkgdir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) 416 if pkgdir not in self.pythonpaths: 417 self.pythonpaths.append(pkgdir)

418 - def make_solver_keywords(self):

419 kw = super(CeseCase, self).make_solver_keywords() 420 kw['debug'] = self.solver.debug_cese 421 for key in 'diffname', 'tauname', 'omeganame': 422 val = self.solver.get(key) 423 if val != None: kw[key] = val 424 kw['alpha'] = int(self.solver.alpha) 425 for key in ('taylor', 'cnbfac', 'sftfac', 426 'taumin', 'taumax', 'tauscale', 'omegamin', 'omegascale', 427 'mqmin', 'mqscale',): 428 val = self.solver.get(key) 429 if val != None: kw[key] = float(val) 430 return kw

431

432 ################################################################################ 433 # Boundary conditions. 434 ################################################################################ 435 436 -class CeseBC(BC):

437 """ 438 Base class for all BC types for CESE method, except periodic BC. 439 440 @cvar _ghostgeom_: selector for the ghost geometry caculator. 441 @ctype _ghostgeom_: str 442 """ 443 from solvcon.dependency import getcdll 444 __clib_ceseb = { 445 2: getcdll('ceseb2d'), 446 3: getcdll('ceseb3d'), 447 } 448 del getcdll 449 @property

450 - def _clib_ceseb(self):

451 return self.__clib_ceseb[self.svr.ndim]

452 453 _ghostgeom_ = None

454 - def init(self, **kw):

455 from ctypes import byref, c_int 456 super(CeseBC, self).init(**kw) 457 getattr(self._clib_ceseb, 'ghostgeom_'+self._ghostgeom_)( 458 byref(self.svr.exd), 459 c_int(self.facn.shape[0]), 460 self.facn.ctypes._as_parameter_, 461 )

462 - def soln(self):

463 """ 464 Update ghost cells after marchsol. 465 """ 466 pass

467 - def dsoln(self):

468 """ 469 Update ghost cells after marchdsol. 470 """ 471 pass

472

473 -class CeseNonrefl(CeseBC):

474 """ 475 General non-reflective boundary condition (NRBC). 476 """ 477 _ghostgeom_ = 'mirror'

478 - def soln(self):

479 from ctypes import byref, c_int 480 svr = self.svr 481 self._clib_ceseb.bound_nonrefl_soln( 482 byref(svr.exd), 483 c_int(self.facn.shape[0]), 484 self.facn.ctypes._as_parameter_, 485 )

486 - def dsoln(self):

487 from ctypes import byref, c_int 488 svr = self.svr 489 self._clib_ceseb.bound_nonrefl_dsoln( 490 byref(svr.exd), 491 c_int(self.facn.shape[0]), 492 self.facn.ctypes._as_parameter_, 493 )

494

495 -class CesePeriodic(periodic):

496 """ 497 General periodic boundary condition for sequential runs. 498 """

499 - def init(self, **kw):

500 svr = self.svr 501 ngstcell = svr.ngstcell 502 ngstface = svr.ngstface 503 facn = self.facn 504 slctm = self.rclp[:,0] + ngstcell 505 slctr = self.rclp[:,1] + ngstcell 506 # move coordinates. 507 shf = svr.cecnd[slctr,0,:] - svr.fccnd[facn[:,2]+ngstface,:] 508 svr.cecnd[slctm,0,:] = svr.fccnd[facn[:,0]+ngstface,:] + shf

509 - def soln(self):

510 svr = self.svr 511 slctm = self.rclp[:,0] + svr.ngstcell 512 slctr = self.rclp[:,1] + svr.ngstcell 513 svr.soln[slctm,:] = svr.soln[slctr,:]

514 - def dsoln(self):

515 svr = self.svr 516 slctm = self.rclp[:,0] + svr.ngstcell 517 slctr = self.rclp[:,1] + svr.ngstcell 518 svr.dsoln[slctm,:,:] = svr.dsoln[slctr,:,:]

519

520 ################################################################################ 521 # Anchors. 522 ################################################################################ 523 524 -class ConvergeAnchor(Anchor):

525 - def __init__(self, svr, **kw):

526 self.norm = {} 527 super(ConvergeAnchor, self).__init__(svr, **kw)

528 - def preloop(self):

529 from numpy import empty 530 svr = self.svr 531 der = svr.der 532 der['diff'] = empty((svr.ngstcell+svr.ncell, svr.neq), 533 dtype=svr.fpdtype)

534 - def postfull(self):

535 from ctypes import c_int, c_double 536 svr = self.svr 537 diff = svr.der['diff'] 538 svr._tcall(svr._clib_cese.calc_norm_diff, -svr.ngstcell, svr.ncell, 539 diff[svr.ngstcell:].ctypes._as_parameter_, 540 ) 541 # Linf norm. 542 Linf = [] 543 Linf.extend(diff.max(axis=0)) 544 self.norm['Linf'] = Linf 545 # L1 norm. 546 svr._clib_cese.calc_norm_L1.restype = c_double 547 L1 = [] 548 for ieq in range(svr.neq): 549 vals = svr._tcall(svr._clib_cese.calc_norm_L1, 0, svr.ncell, 550 diff[svr.ngstcell:].ctypes._as_parameter_, 551 c_int(ieq), 552 ) 553 L1.append(sum(vals)) 554 self.norm['L1'] = L1 555 # L2 norm. 556 svr._clib_cese.calc_norm_L2.restype = c_double 557 L2 = [] 558 for ieq in range(svr.neq): 559 vals = svr._tcall(svr._clib_cese.calc_norm_L2, 0, svr.ncell, 560 diff[svr.ngstcell:].ctypes._as_parameter_, 561 c_int(ieq), 562 ) 563 L2.append(sum(vals)) 564 self.norm['L2'] = L2

565

566 -class Probe(object):

567 """ 568 Represent a point in the mesh. 569 """

570 - def __init__(self, *args, **kw):

571 from numpy import array 572 self.speclst = kw.pop('speclst') 573 self.name = kw.pop('name', None) 574 self.crd = array(args, dtype='float64') 575 self.pcl = -1 576 self.vals = list()

577 - def __str__(self):

578 crds = ','.join(['%g'%val for val in self.crd]) 579 return 'Pt/%s#%d(%s)%d' % (self.name, self.pcl, crds, len(self.vals))

580 - def locate_cell(self, svr):

581 idx = svr.locate_point(*self.crd) 582 self.pcl = idx[0]

583 - def __call__(self, svr, time):

584 ngstcell = svr.ngstcell 585 vlist = [time] 586 for spec in self.speclst: 587 arr = None 588 if isinstance(spec, str): 589 arr = svr.der[spec] 590 elif isinstance(spec, int): 591 if spec >= 0 and spec < svr.neq: 592 arr = svr.soln[:,spec] 593 elif spec < 0 and -1-spec < svr.neq: 594 spec = -1-spec 595 arr = svr.sol[:,spec] 596 if arr == None: 597 raise IndexError, 'spec %s incorrect'%str(spec) 598 vlist.append(arr[ngstcell+self.pcl]) 599 self.vals.append(vlist)

600

601 -class ProbeAnchor(Anchor):

602 """ 603 Anchor for probe. 604 """

605 - def __init__(self, svr, **kw):

606 speclst = kw.pop('speclst') 607 self.points = list() 608 for data in kw.pop('coords'): 609 pkw = {'speclst': speclst, 'name': data[0]} 610 self.points.append(Probe(*data[1:], **pkw)) 611 super(ProbeAnchor, self).__init__(svr, **kw)

612 - def preloop(self):

613 for point in self.points: point.locate_cell(self.svr) 614 for point in self.points: point(self.svr, self.svr.time)

615 - def postfull(self):

616 for point in self.points: point(self.svr, self.svr.time)

617

618 ################################################################################ 619 # Hooks. 620 ################################################################################ 621 622 -class CflHook(Hook):

623 """ 624 Make sure is CFL number is bounded and print averaged CFL number over time. 625 626 @ivar cflmin: CFL number should be greater than or equal to the value. 627 @itype cflmin: float 628 @ivar cflmax: CFL number should be less than the value. 629 @itype cflmax: float 630 @ivar fullstop: flag to stop when CFL is out of bound. Default True. 631 @itype fullstop: bool 632 """

633 - def __init__(self, cse, **kw):

634 self.cflmin = kw.pop('cflmin', 0.0) 635 self.cflmax = kw.pop('cflmax', 1.0) 636 self.fullstop = kw.pop('fullstop', True) 637 self.aCFL = 0.0 638 self.mCFL = 0.0 639 self.hnCFL = 1.0 640 self.hxCFL = 0.0 641 self.aadj = 0 642 self.haadj = 0 643 super(CflHook, self).__init__(cse, **kw)

644 - def _notify(self, msg):

645 from warnings import warn 646 if self.fullstop: 647 raise RuntimeError, msg 648 else: 649 warn(msg)

650 - def postmarch(self):

651 from numpy import isnan 652 info = self.info 653 steps_stride = self.cse.execution.steps_stride 654 istep = self.cse.execution.step_current 655 mr = self.cse.execution.marchret 656 isp = self.cse.is_parallel 657 psteps = self.psteps 658 # collect CFL. 659 nCFL = max([m['cfl'][0] for m in mr]) if isp else mr['cfl'][0] 660 xCFL = max([m['cfl'][1] for m in mr]) if isp else mr['cfl'][1] 661 nadj = sum([m['cfl'][2] for m in mr]) if isp else mr['cfl'][2] 662 aadj = sum([m['cfl'][3] for m in mr]) if isp else mr['cfl'][2] 663 hnCFL = min([nCFL, self.hnCFL]) 664 self.hnCFL = hnCFL if not isnan(hnCFL) else self.hnCFL 665 hxCFL = max([xCFL, self.hxCFL]) 666 self.hxCFL = hxCFL if not isnan(hxCFL) else self.hxCFL 667 self.aCFL += xCFL*steps_stride 668 self.mCFL = self.aCFL/(istep+steps_stride) 669 self.aadj += aadj 670 self.haadj += aadj 671 # check. 672 if self.cflmin != None and nCFL < self.cflmin: 673 self._notify("CFL = %g < %g after step: %d" % ( 674 nCFL, self.cflmin, istep)) 675 if self.cflmax != None and xCFL >= self.cflmax: 676 self._notify("CFL = %g >= %g after step: %d" % ( 677 xCFL, self.cflmax, istep)) 678 # output information. 679 if istep > 0 and istep%psteps == 0: 680 info("CFL = %.2f/%.2f - %.2f/%.2f adjusted: %d/%d/%d\n" % ( 681 nCFL, xCFL, self.hnCFL, self.hxCFL, nadj, 682 self.aadj, self.haadj)) 683 self.aadj = 0

684

685 - def postloop(self):

686 self.info("Averaged maximum CFL = %g.\n" % self.mCFL)

687

688 -class ConvergeHook(BlockHook):

689 - def __init__(self, cse, **kw):

690 self.name = kw.pop('name', 'converge') 691 self.keys = kw.pop('keys', None) 692 self.eqs = kw.pop('eqs', None) 693 csteps = kw.pop('csteps', None) 694 super(ConvergeHook, self).__init__(cse, **kw) 695 self.csteps = self.psteps if csteps == None else cstep 696 self.ankkw = kw

697 - def drop_anchor(self, svr):

698 ankkw = self.ankkw.copy() 699 ankkw['name'] = self.name 700 self._deliver_anchor(svr, ConvergeAnchor, ankkw)

701 - def _collect(self):

702 from numpy import sqrt 703 cse = self.cse 704 neq = cse.execution.neq 705 if cse.is_parallel: 706 dom = cse.solver.domainobj 707 dealer = cse.solver.dealer 708 allnorm = list() 709 for iblk in range(dom.nblk): 710 dealer[iblk].cmd.pullank(self.name, 'norm', with_worker=True) 711 allnorm.append(dealer[iblk].recv()) 712 norm = {'Linf': [0.0]*neq, 'L1': [0.0]*neq, 'L2': [0.0]*neq} 713 for ieq in range(neq): 714 norm['Linf'][ieq] = max([nm['Linf'][ieq] for nm in allnorm]) 715 norm['L1'][ieq] = sum([nm['L1'][ieq] for nm in allnorm]) 716 norm['L2'][ieq] = sum([nm['L2'][ieq] for nm in allnorm]) 717 else: 718 svr = self.cse.solver.solverobj 719 norm = svr.runanchors[self.name].norm 720 for ieq in range(neq): 721 norm['L2'][ieq] = sqrt(norm['L2'][ieq]) 722 self.norm = norm

723 - def postmarch(self):

724 info = self.info 725 cse = self.cse 726 istep = cse.execution.step_current 727 csteps = self.csteps 728 psteps = self.psteps 729 if istep > 0 and istep%csteps == 0: 730 self._collect() 731 if istep > 0 and istep%psteps == 0: 732 norm = self.norm 733 keys = self.keys if self.keys != None else self.norm.keys() 734 keys.sort() 735 eqs = self.eqs if self.eqs != None else range(cse.execution.neq) 736 for key in keys: 737 info("Converge/%-4s [ %s ]:\n %s\n" % (key, 738 ', '.join(['%d'%ieq for ieq in eqs]), 739 ' '.join(['%.4e'%norm[key][ieq] for ieq in eqs]), 740 ))

741

742 -class VerifyHook(BlockHook):

743 - def __init__(self, cse, **kw):

744 import os 745 self.arrnames = kw.pop('arrnames') 746 self.goldenext = kw.pop('goldenext', 'golden.npy') 747 self.goldendir = os.path.abspath(kw.pop('goldendir', 748 os.path.join(cse.io.rootdir, 'golden') 749 )) 750 super(VerifyHook, self).__init__(cse, **kw)

751 - def _verify(self):

752 import os 753 from numpy import save, load, abs 754 for arrname in self.arrnames: 755 arr = self.cse.execution.var[arrname] 756 goldenfn = '.'.join([self.cse.io.basefn, arrname, self.goldenext]) 757 goldenpath = os.path.join(self.goldendir, goldenfn) 758 if os.path.exists(goldenpath): 759 golden = load(goldenpath) 760 diffsum = abs(arr - golden).sum() 761 self.cse.execution.verified_norm = diffsum 762 self.info('Verify %s: %e' % (arrname, diffsum)) 763 if diffsum == 0: 764 self.info(': good.\n') 765 else: 766 self.info(': BAD!\n') 767 else: 768 goldenpath = os.path.join(self.cse.io.basedir, goldenfn) 769 self.info('Save %s to %s .\n' % (arrname, goldenpath)) 770 save(goldenpath, arr)

771 - def postloop(self):

772 self._verify()

773

774 -class ProbeHook(BlockHook):

775 """ 776 Point probe. 777 """

778 - def __init__(self, cse, **kw):

779 self.name = kw.pop('name', 'ppank') 780 super(ProbeHook, self).__init__(cse, **kw) 781 self.ankkw = kw 782 self.points = None

783 - def drop_anchor(self, svr):

784 ankkw = self.ankkw.copy() 785 ankkw['name'] = self.name 786 self._deliver_anchor(svr, ProbeAnchor, ankkw)

787 - def _collect(self):

788 cse = self.cse 789 if cse.is_parallel: 790 dom = cse.solver.domainobj 791 dealer = cse.solver.dealer 792 allpoints = list() 793 for iblk in range(dom.nblk): 794 dealer[iblk].cmd.pullank(self.name, 'points', with_worker=True) 795 allpoints.append(dealer[iblk].recv()) 796 npt = len(allpoints[0]) 797 points = [None]*npt 798 for rpoints in allpoints: 799 ipt = 0 800 while ipt < npt: 801 if points[ipt] == None and rpoints[ipt].pcl >=0: 802 points[ipt] = rpoints[ipt] 803 ipt += 1 804 else: 805 svr = self.cse.solver.solverobj 806 points = [pt for pt in svr.runanchors[self.name].points 807 if pt.pcl >= 0] 808 self.points = points

809 - def postmarch(self):

810 psteps = self.psteps 811 istep = self.cse.execution.step_current 812 if istep%psteps != 0: return False 813 self._collect() 814 return True

815 - def postloop(self):

816 import os 817 from numpy import array, save 818 for point in self.points: 819 ptfn = '%s_pt_%s_%s.npy' % ( 820 self.cse.io.basefn, self.name, point.name) 821 ptfn = os.path.join(self.cse.io.basedir, ptfn) 822 save(ptfn, array(point.vals, dtype='float64'))

823

Source Code for Module solvcon.kerpak.cese