solvcon.kerpak.lincese

1 # -*- coding: UTF-8 -*- 2 # 3 # Copyright (C) 2010-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 CESE solver specialized for general linear equations. 21 """ 22 23 from solvcon.gendata import SingleAssignDict, AttributeDict 24 from solvcon.anchor import Anchor 25 from solvcon.hook import BlockHook 26 from solvcon.kerpak.cese import CeseSolver, CeseCase

27 28 ############################################################################### 29 # Solver. 30 ############################################################################### 31 32 -class LinceseSolver(CeseSolver):

33 """ 34 Basic linear CESE solver. 35 36 @ivar cfldt: the time_increment for CFL calculation at boundcond. 37 @itype cfldt: float 38 @ivar cflmax: the maximum CFL number. 39 @itype cflmax: float 40 """ 41 from solvcon.dependency import getcdll 42 __clib_lincese = { 43 2: getcdll('lincese2d'), 44 3: getcdll('lincese3d'), 45 } 46 del getcdll 47 @property

48 - def _clib_lincese(self):

49 return self.__clib_lincese[self.ndim]

50 @property

51 - def _jacofunc_(self):

52 return self._clib_lincese.calc_jaco

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

54 self.cfldt = kw.pop('cfldt', None) 55 self.cflmax = 0.0 56 super(LinceseSolver, self).__init__(*args, **kw)

57 - def make_grpda(self):

58 raise NotImplementedError

59 - def provide(self):

60 from ctypes import byref, c_int 61 # fill group data array. 62 self.make_grpda() 63 # pre-calculate CFL. 64 self._set_time(self.time, self.cfldt) 65 self._clib_lincese.calc_cfl( 66 byref(self.exd), c_int(0), c_int(self.ncell)) 67 self.cflmax = self.cfl.max() 68 # super method. 69 super(LinceseSolver, self).provide()

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

71 self.marchret.setdefault('cfl', [0.0, 0.0, 0, 0]) 72 self.marchret['cfl'][0] = self.cflmax 73 self.marchret['cfl'][1] = self.cflmax 74 self.marchret['cfl'][2] = 0 75 self.marchret['cfl'][3] = 0 76 return self.marchret

77

78 ############################################################################### 79 # Case. 80 ############################################################################### 81 82 -class LinceseCase(CeseCase):

83 """ 84 Basic case with linear CESE method. 85 """ 86 from solvcon.domain import Domain 87 defdict = { 88 'solver.solvertype': LinceseSolver, 89 'solver.domaintype': Domain, 90 'solver.alpha': 0, 91 'solver.cfldt': None, 92 } 93 del Domain

94 - def make_solver_keywords(self):

95 kw = super(LinceseCase, self).make_solver_keywords() 96 # setup delta t for CFL calculation. 97 cfldt = self.solver.cfldt 98 cfldt = self.execution.time_increment if cfldt is None else cfldt 99 kw['cfldt'] = cfldt 100 return kw

101

102 ############################################################################### 103 # Plane wave solution and initializer. 104 ############################################################################### 105 106 -class PlaneWaveSolution(object):

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

108 from numpy import sqrt 109 wvec = kw['wvec'] 110 ctr = kw['ctr'] 111 amp = kw['amp'] 112 assert len(wvec) == len(ctr) 113 # calculate eigenvalues and eigenvectors. 114 evl, evc = self._calc_eigen(**kw) 115 # store data to self. 116 self.amp = evc * (amp / sqrt((evc**2).sum())) 117 self.ctr = ctr 118 self.wvec = wvec 119 self.afreq = evl * sqrt((wvec**2).sum()) 120 self.wsp = evl

121 - def _calc_eigen(self, *args, **kw):

122 """ 123 Calculate eigenvalues and eigenvectors. 124 125 @return: eigenvalues and eigenvectors. 126 @rtype: tuple 127 """ 128 raise NotImplementedError

129 - def __call__(self, svr, asol, adsol):

130 from ctypes import byref, c_double 131 svr._clib_lincese.calc_planewave( 132 byref(svr.exd), 133 asol.ctypes._as_parameter_, 134 adsol.ctypes._as_parameter_, 135 self.amp.ctypes._as_parameter_, 136 self.ctr.ctypes._as_parameter_, 137 self.wvec.ctypes._as_parameter_, 138 c_double(self.afreq), 139 )

140

141 -class PlaneWaveAnchor(Anchor):

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

143 self.planewaves = kw.pop('planewaves') 144 super(PlaneWaveAnchor, self).__init__(svr, **kw)

145 - def _calculate(self, asol):

146 for pw in self.planewaves: 147 pw(self.svr, asol, self.adsol)

148 - def provide(self):

149 from numpy import empty 150 ngstcell = self.svr.ngstcell 151 nacell = self.svr.ncell + ngstcell 152 # plane wave solution. 153 asol = self.svr.der['analytical'] = empty( 154 (nacell, self.svr.neq), dtype=self.svr.fpdtype) 155 adsol = self.adsol = empty( 156 (nacell, self.svr.neq, self.svr.ndim), 157 dtype=self.svr.fpdtype) 158 asol.fill(0.0) 159 self._calculate(asol) 160 self.svr.soln[ngstcell:,:] = asol[ngstcell:,:] 161 self.svr.dsoln[ngstcell:,:,:] = adsol[ngstcell:,:,:] 162 # difference. 163 diff = self.svr.der['difference'] = empty( 164 (nacell, self.svr.neq), dtype=self.svr.fpdtype) 165 diff[ngstcell:,:] = self.svr.soln[ngstcell:,:] - asol[ngstcell:,:]

166 - def postfull(self):

167 ngstcell = self.svr.ngstcell 168 # plane wave solution. 169 asol = self.svr.der['analytical'] 170 asol.fill(0.0) 171 self._calculate(asol) 172 # difference. 173 diff = self.svr.der['difference'] 174 diff[ngstcell:,:] = self.svr.soln[ngstcell:,:] - asol[ngstcell:,:]

175

176 -class PlaneWaveHook(BlockHook):

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

178 self.planewaves = kw.pop('planewaves') 179 self.norm = dict() 180 super(PlaneWaveHook, self).__init__(svr, **kw)

181 - def drop_anchor(self, svr):

182 svr.runanchors.append( 183 PlaneWaveAnchor(svr, planewaves=self.planewaves) 184 )

185 - def _calculate(self):

186 from numpy import empty, sqrt, abs 187 neq = self.cse.execution.neq 188 var = self.cse.execution.var 189 asol = self._collect_interior( 190 'analytical', inder=True, consider_ghost=True) 191 diff = self._collect_interior( 192 'difference', inder=True, consider_ghost=True) 193 norm_Linf = empty(neq, dtype='float64') 194 norm_L2 = empty(neq, dtype='float64') 195 clvol = self.blk.clvol 196 for it in range(neq): 197 norm_Linf[it] = abs(diff[:,it]).max() 198 norm_L2[it] = sqrt((diff[:,it]**2*clvol).sum()) 199 self.norm['Linf'] = norm_Linf 200 self.norm['L2'] = norm_L2

201 - def preloop(self):

202 from numpy import pi 203 self.postmarch() 204 for ipw in range(len(self.planewaves)): 205 pw = self.planewaves[ipw] 206 self.info("planewave[%d]:\n" % ipw) 207 self.info(" c = %g, omega = %g, T = %.15e\n" % ( 208 pw.wsp, pw.afreq, 2*pi/pw.afreq))

209 - def postmarch(self):

210 psteps = self.psteps 211 istep = self.cse.execution.step_current 212 if istep%psteps == 0: 213 self._calculate()

214 - def postloop(self):

215 import os 216 from cPickle import dump 217 fname = '%s_norm.pickle' % self.cse.io.basefn 218 fname = os.path.join(self.cse.io.basedir, fname) 219 dump(self.norm, open(fname, 'wb'), -1) 220 self.info('Linf norm in velocity:\n') 221 self.info(' %e, %e, %e\n' % tuple(self.norm['Linf'][:3])) 222 self.info('L2 norm in velocity:\n') 223 self.info(' %e, %e, %e\n' % tuple(self.norm['L2'][:3]))

224

Source Code for Module solvcon.kerpak.lincese