solvcon.kerpak.euler

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 Euler equations solver using the CESE method. 21 """ 22 23 from .cese import CeseSolver 24 from .cese import CeseCase 25 from .cese import CeseBC 26 from solvcon.anchor import Anchor 27 from solvcon.hook import BlockHook

28 29 ############################################################################### 30 # Solver. 31 ############################################################################### 32 33 -class EulerSolver(CeseSolver):

34 """ 35 Inviscid aerodynamic solver for the Euler equations. 36 """

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

38 self.cflname = kw.pop('cflname', 'adj') 39 kw['nsca'] = 1 40 super(EulerSolver, self).__init__(blk, *args, **kw) 41 self.cflc = self.cfl.copy() # FIXME: obselete?

42 from solvcon.dependency import getcdll 43 __clib_euler = { 44 2: getcdll('euler2d'), 45 3: getcdll('euler3d'), 46 } 47 del getcdll 48 @property

49 - def _clib_euler(self):

50 return self.__clib_euler[self.ndim]

51 _gdlen_ = 0 52 @property

53 - def _jacofunc_(self):

54 return self._clib_euler.calc_jaco

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

56 func = getattr(self._clib_euler, 'calc_cfl_'+self.cflname) 57 self._tcall(func, 0, self.ncell) 58 mincfl = self.ocfl.min() 59 maxcfl = self.ocfl.max() 60 nadj = (self.cfl==1).sum() 61 self.marchret.setdefault('cfl', [0.0, 0.0, 0, 0]) 62 self.marchret['cfl'][0] = mincfl 63 self.marchret['cfl'][1] = maxcfl 64 self.marchret['cfl'][2] = nadj 65 self.marchret['cfl'][3] += nadj 66 return self.marchret

67

68 ############################################################################### 69 # Case. 70 ############################################################################### 71 72 -class EulerCase(CeseCase):

73 """ 74 Inviscid aerodynamic case for the Euler equations. 75 """ 76 from solvcon.domain import Domain 77 defdict = { 78 'solver.solvertype': EulerSolver, 79 'solver.domaintype': Domain, 80 'solver.cflname': 'adj', 81 } 82 del Domain

83 - def make_solver_keywords(self):

84 kw = super(EulerCase, self).make_solver_keywords() 85 kw['cflname'] = self.solver.cflname 86 return kw

87 - def load_block(self):

88 loaded = super(EulerCase, self).load_block() 89 if hasattr(loaded, 'ndim'): 90 ndim = loaded.ndim 91 else: 92 ndim = loaded.blk.ndim 93 self.execution.neq = ndim+2 94 return loaded

95

96 ############################################################################### 97 # Boundary conditions. 98 ############################################################################### 99 100 -class EulerBC(CeseBC):

101 """ 102 Basic BC class for the Euler equations. 103 """ 104 from solvcon.dependency import getcdll 105 __clib_eulerb = { 106 2: getcdll('eulerb2d'), 107 3: getcdll('eulerb3d'), 108 } 109 del getcdll 110 @property

111 - def _clib_eulerb(self):

112 return self.__clib_eulerb[self.svr.ndim]

113

114 -class EulerWall(EulerBC):

115 _ghostgeom_ = 'mirror'

116 - def soln(self):

117 from ctypes import byref, c_int 118 svr = self.svr 119 self._clib_eulerb.bound_wall_soln( 120 byref(svr.exd), 121 c_int(self.facn.shape[0]), 122 self.facn.ctypes._as_parameter_, 123 )

124 - def dsoln(self):

125 from ctypes import byref, c_int 126 svr = self.svr 127 self._clib_eulerb.bound_wall_dsoln( 128 byref(svr.exd), 129 c_int(self.facn.shape[0]), 130 self.facn.ctypes._as_parameter_, 131 )

132 -class EulerNonslipWall(EulerBC):

133 _ghostgeom_ = 'mirror'

134 - def soln(self):

135 from ctypes import byref, c_int 136 svr = self.svr 137 self._clib_eulerb.bound_nonslipwall_soln( 138 byref(svr.exd), 139 c_int(self.facn.shape[0]), 140 self.facn.ctypes._as_parameter_, 141 )

142 - def dsoln(self):

143 from ctypes import byref, c_int 144 svr = self.svr 145 self._clib_eulerb.bound_nonslipwall_dsoln( 146 byref(svr.exd), 147 c_int(self.facn.shape[0]), 148 self.facn.ctypes._as_parameter_, 149 )

150

151 -class EulerInlet(EulerBC):

152 vnames = ['rho', 'v1', 'v2', 'v3', 'p', 'gamma'] 153 vdefaults = { 154 'rho': 1.0, 'p': 1.0, 'gamma': 1.4, 'v1': 0.0, 'v2': 0.0, 'v3': 0.0, 155 } 156 _ghostgeom_ = 'mirror'

157 - def soln(self):

158 from ctypes import byref, c_int 159 svr = self.svr 160 self._clib_eulerb.bound_inlet_soln( 161 byref(svr.exd), 162 c_int(self.facn.shape[0]), 163 self.facn.ctypes._as_parameter_, 164 c_int(self.value.shape[1]), 165 self.value.ctypes._as_parameter_, 166 )

167 - def dsoln(self):

168 from ctypes import byref, c_int 169 svr = self.svr 170 self._clib_eulerb.bound_inlet_dsoln( 171 byref(svr.exd), 172 c_int(self.facn.shape[0]), 173 self.facn.ctypes._as_parameter_, 174 )

175

176 -class EulerOutlet(EulerBC):

177 _ghostgeom_ = 'translate'

178 - def soln(self):

179 from ctypes import byref, c_int 180 svr = self.svr 181 self._clib_eulerb.bound_outlet_soln( 182 byref(svr.exd), 183 c_int(self.facn.shape[0]), 184 self.facn.ctypes._as_parameter_, 185 )

186 - def dsoln(self):

187 from ctypes import byref, c_int 188 svr = self.svr 189 self._clib_eulerb.bound_outlet_dsoln( 190 byref(svr.exd), 191 c_int(self.facn.shape[0]), 192 self.facn.ctypes._as_parameter_, 193 )

194

195 ############################################################################### 196 # Anchors. 197 ############################################################################### 198 199 -class EulerIAnchor(Anchor):

200 """ 201 Basic initializing anchor class for all Euler problems. 202 """

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

204 assert isinstance(svr, EulerSolver) 205 self.gamma = float(kw.pop('gamma')) 206 super(EulerIAnchor, self).__init__(svr, **kw)

207 - def provide(self):

208 from solvcon.solver import ALMOST_ZERO 209 svr = self.svr 210 svr.amsca.fill(self.gamma) 211 svr.sol.fill(ALMOST_ZERO) 212 svr.soln.fill(ALMOST_ZERO) 213 svr.dsol.fill(ALMOST_ZERO) 214 svr.dsoln.fill(ALMOST_ZERO)

215

216 -class UniformIAnchor(EulerIAnchor):

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

218 self.rho = float(kw.pop('rho')) 219 self.v1 = float(kw.pop('v1')) 220 self.v2 = float(kw.pop('v2')) 221 self.v3 = float(kw.pop('v3')) 222 self.p = float(kw.pop('p')) 223 super(UniformIAnchor, self).__init__(svr, **kw)

224 - def provide(self):

225 super(UniformIAnchor, self).provide() 226 gamma = self.gamma 227 svr = self.svr 228 svr.soln[:,0].fill(self.rho) 229 svr.soln[:,1].fill(self.rho*self.v1) 230 svr.soln[:,2].fill(self.rho*self.v2) 231 vs = self.v1**2 + self.v2**2 232 if svr.ndim == 3: 233 vs += self.v3**2 234 svr.soln[:,3].fill(self.rho*self.v3) 235 svr.soln[:,svr.ndim+1].fill(self.rho*vs/2 + self.p/(gamma-1)) 236 svr.sol[:] = svr.soln[:]

237

238 -class EulerOAnchor(Anchor):

239 _varlist_ = ['v', 'rho', 'p', 'T', 'ke', 'a', 'M', 'sch']

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

241 self.schk = kw.pop('schk', 1.0) 242 self.schk0 = kw.pop('schk0', 0.0) 243 self.schk1 = kw.pop('schk1', 1.0) 244 super(EulerOAnchor, self).__init__(svr, **kw)

245 - def provide(self):

246 from numpy import empty 247 svr = self.svr 248 der = svr.der 249 nelm = svr.ngstcell + svr.ncell 250 der['v'] = empty((nelm, svr.ndim), dtype=svr.fpdtype) 251 der['rho'] = empty(nelm, dtype=svr.fpdtype) 252 der['p'] = empty(nelm, dtype=svr.fpdtype) 253 der['T'] = empty(nelm, dtype=svr.fpdtype) 254 der['ke'] = empty(nelm, dtype=svr.fpdtype) 255 der['a'] = empty(nelm, dtype=svr.fpdtype) 256 der['M'] = empty(nelm, dtype=svr.fpdtype) 257 der['sch'] = empty(nelm, dtype=svr.fpdtype) 258 self._calculate_physics() 259 self._calculate_schlieren()

260 - def _calculate_physics(self):

261 svr = self.svr 262 der = svr.der 263 svr._tcall(svr._clib_euler.calc_physics, -svr.ngstcell, svr.ncell, 264 der['v'].ctypes._as_parameter_, 265 der['rho'].ctypes._as_parameter_, 266 der['p'].ctypes._as_parameter_, 267 der['T'].ctypes._as_parameter_, 268 der['ke'].ctypes._as_parameter_, 269 der['a'].ctypes._as_parameter_, 270 der['M'].ctypes._as_parameter_, 271 )

272 - def _calculate_schlieren(self):

273 from ctypes import c_double 274 svr = self.svr 275 sch = svr.der['sch'] 276 svr._tcall(svr._clib_euler.calc_schlieren_rhog, 277 -svr.ngstcell, svr.ncell, sch.ctypes._as_parameter_) 278 rhogmax = sch[svr.ngstcell:].max() 279 svr._tcall(svr._clib_euler.calc_schlieren_sch, 280 -svr.ngstcell, svr.ncell, 281 c_double(self.schk), c_double(self.schk0), c_double(self.schk1), 282 c_double(rhogmax), sch.ctypes._as_parameter_, 283 )

284 - def postfull(self):

285 self._calculate_physics() 286 self._calculate_schlieren()

287

288 ############################################################################### 289 # Hooks. 290 ############################################################################### 291 292 -class ConservationHook(BlockHook):

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

294 self.records = list() 295 self.recordpath = None 296 super(ConservationHook, self).__init__(cse, **kw)

297 - def preloop(self):

298 import os 299 recordfn = '%s_conservation.npy' % self.cse.io.basefn 300 self.recordpath = os.path.join(self.cse.io.basedir, recordfn)

301 - def _calculate(self):

302 clvol = self.blk.clvol 303 rho = self.cse.execution.var['rho'] 304 etot = self.cse.execution.var['soln'][:,-1] 305 utot = etot - self.cse.execution.var['ke'] 306 crho = (clvol*rho).sum() 307 cetot = (clvol*etot).sum() 308 cutot = (clvol*utot).sum() 309 self.records.append([self.cse.execution.time, 310 crho, cetot, cutot])

311 - def postmarch(self):

312 from numpy import array, save 313 self._calculate() 314 istep = self.cse.execution.step_current 315 if istep%self.psteps == 0: 316 info = self.info 317 d0 = array(self.records[0][1:], dtype='float64') 318 dn = array(self.records[-1][1:], dtype='float64') 319 diff = (dn-d0)/d0 * 100 320 info('density conservation difference: %7.3f%%\n' % diff[0]) 321 info('energy conservation difference: %7.3f%%\n' % diff[1]) 322 info('thermal conservation difference: %7.3f%%\n' % diff[2]) 323 save(self.recordpath, self.records)

324 - def postloop(self):

325 from numpy import array 326 info = self.info 327 self.records = array(self.records, dtype='float64') 328 info('mass conservation: %g %g\n' % ( 329 self.records[0,1], self.records[-1,1])) 330 info('energy conservation: %g %g\n' % ( 331 self.records[0,2], self.records[-1,2])) 332 info('thermal conservation: %g %g\n' % ( 333 self.records[0,3], self.records[-1,3])) 334 info('Conservation record saved to %s .\n' % self.recordpath)

335

Source Code for Module solvcon.kerpak.euler