solvcon.kerpak.lincuse

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 with CUDA support. 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.cuse import CUDA_RAISE_ON_FAIL, CuseSolver, CuseCase

27 28 ############################################################################### 29 # Solver. 30 ############################################################################### 31 32 -class LincuseSolver(CuseSolver):

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_lincuse_c = { 43 2: getcdll('lincuse2d_c'), 44 3: getcdll('lincuse3d_c'), 45 } 46 __clib_lincuse_cu = { 47 2: getcdll('lincuse2d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 48 3: getcdll('lincuse3d_cu', raise_on_fail=CUDA_RAISE_ON_FAIL), 49 } 50 del getcdll 51 @property

52 - def _clib_lincuse_c(self):

53 return self.__clib_lincuse_c[self.ndim]

54 @property

55 - def _clib_lincuse_cu(self):

56 return self.__clib_lincuse_cu[self.ndim]

57 @property

58 - def _jacofunc_(self):

59 return self._clib_lincuse_c.calc_jaco

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

61 self.cfldt = kw.pop('cfldt', None) 62 self.cflmax = 0.0 63 super(LincuseSolver, self).__init__(*args, **kw)

64 - def make_grpda(self):

65 raise NotImplementedError

66 - def provide(self):

67 from ctypes import byref, c_int 68 # fill group data array. 69 self.make_grpda() 70 # pre-calculate CFL. 71 self._set_time(self.time, self.cfldt) 72 self._clib_lincuse_c.calc_cfl( 73 byref(self.exd), c_int(0), c_int(self.ncell)) 74 self.ocfl[:] = self.cfl[:] 75 # super method. 76 super(LincuseSolver, self).provide()

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

78 pass

79

80 ############################################################################### 81 # Case. 82 ############################################################################### 83 84 -class LincuseCase(CuseCase):

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

96 - def make_solver_keywords(self):

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

103

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

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

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

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

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

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

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

142

143 -class PlaneWaveAnchor(Anchor):

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

145 self.planewaves = kw.pop('planewaves') 146 super(PlaneWaveAnchor, self).__init__(svr, **kw)

147 - def _calculate(self, asol):

148 for pw in self.planewaves: 149 pw(self.svr, asol, self.adsol)

150 - def provide(self):

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

168 - def postfull(self):

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

177

178 -class PlaneWaveHook(BlockHook):

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

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

183 - def drop_anchor(self, svr):

184 svr.runanchors.append( 185 PlaneWaveAnchor(svr, planewaves=self.planewaves) 186 )

187 - def _calculate(self):

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

203 - def preloop(self):

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

211 - def postmarch(self):

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

216 - def postloop(self):

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

226

Source Code for Module solvcon.kerpak.lincuse