######################################################################## # # # ######################################################################## import numpy as np #----------------------------------------------------------------------- # # # #----------------------------------------------------------------------- class wave_system: '''This class defines methods and updates for the scalar advection equation. ''' #------------------------------------------------------------------- # public methods #------------------------------------------------------------------- # The values xmin and xmax need to be modified if one is picky about # periodic boundary conditions. def __init__(self,xmin,xmax,Nx,add_ghost_zones=False): self._N = Nx if add_ghost_zones == True: ng = 3 self._dx = float(xmax - xmin)/(Nx-1-ng) self._xmin = xmin - 2.0*self._dx self._xmax = xmax + self._dx else: self._dx = float(xmax - xmin)/(Nx-1) self._xmin = xmin self._xmax = xmax self._u_n = np.zeros(Nx) self._u_np1 = np.zeros(Nx) self._f_n = np.zeros(Nx) self._dt = 0.1 * self._dx self._flux_type = 1 self._time = 0.0 self._bcond = 1 print "dx = " + str(self._dx) print "dt = " + str(self._dt) #------------------------------------------------------------------- def set_u_n(self, ivp, args): self._u_n = ivp(self.x(), args) # self._u_n[self._N-1] = self._u_n[0] #------------------------------------------------------------------- def cal_flux(self): if (self._flux_type == 1): flux_ftbs(self._f_n, self._u_n, self._N) elif (self._flux_type == 2): flux_ftcs(self._f_n, self._u_n, self._N) elif (self._flux_type == 3): cfl = self._dt/self._dx LaxF(self._f_n, self._u_n, cfl, self._N) elif (self._flux_type == 4): cfl = self._dt/self._dx LaxWendroff_advection(self._f_n, self._u_n, cfl, self._N) elif (self._flux_type == 5): cfl = self._dt/self._dx BeamWarming_advection(self._f_n, self._u_n, cfl, self._N) else: print "Unknown flux_type ",self._flux_type #------------------------------------------------------------------- def set_boundary_condition(self, bctype): self._bcond = bctype #------------------------------------------------------------------- def set_cfl(self, cfl): self._dt = cfl * self._dx #------------------------------------------------------------------- def set_flux_type(self, ft): self._flux_type = ft #------------------------------------------------------------------- def x(self): return np.linspace(self._xmin, self._xmax, self._N) #------------------------------------------------------------------- def u_n(self): return self._u_n[:] #------------------------------------------------------------------- def u_np1(self): return self._u_np1[:] #------------------------------------------------------------------- def f_n(self): return self._f_n[:] #------------------------------------------------------------------- def dx(self): return self._dx #------------------------------------------------------------------- def dt(self): return self._dt #------------------------------------------------------------------- def time(self): return self._time #------------------------------------------------------------------- def flux_type(self): return self._flux_type #------------------------------------------------------------------- def step(self, nsteps): N = self._N dx = self._dx dt = self._dt cfl = dt/dx # up1 = self._u_np1 # u = self._u_n for ii in range (0,nsteps): # print "...taking a step " + str(ii) + " cfl = " + str(cfl) self.cal_flux() for i in range(2,N-2): self._u_np1[i] = self._u_n[i] - cfl*(self._f_n[i] - self._f_n[i-1]) # print "...applying boundary condition" # self._u_np1[0] = self._u_n[0] - cfl*(self._u_n[0] - self._u_n[N-1]) # self._u_np1[0] = self._u_np1[2] # self._u_np1[1] = self._u_np1[2] # self._u_np1[N-2] = self._u_np1[N-3] # self._u_np1[N-1] = self._u_np1[N-3] set_bc(self._u_np1, self._bcond) #rotate pointers temp = self._u_n; self._u_n = self._u_np1; self._u_np1 = temp self._time += self._dt print "...time now " + str(self._time) #------------------------------------------------------------------- # def set_bc(f): # if (self._bcond == 1): # N = size(f) # print "applying outflow bc's. N = " + str(N) # self._u_np1[0] = self._u_np1[2] # self._u_np1[1] = self._u_np1[2] # self._u_np1[N-2] = self._u_np1[N-3] # self._u_np1[N-1] = self._u_np1[N-3] # f[0] = f[2] # f[1] = f[2] # f[N-2] = f[N-3] # f[N-1] = f[N-3] #------------------------------------------------------------------- # Flux routines #------------------------------------------------------------------- def flux_ftbs(f, u, n): for i in range(0,n-1): f[i] = u[i] f[n-1] = f[n-2] def flux_ftcs(f, u, n): for i in range(0,n-1): f[i] = 0.5*(u[i+1] + u[i]) f[n-1] = f[n-2] def LaxF(f, u, cfl, n): for i in range(0,n-1): f[i] = 0.5*(u[i+1] + u[i]) - 0.5*(u[i+1] - u[i])/cfl f[n-1] = f[n-2] def LaxWendroff_advection(f, u, cfl, n): for i in range(0,n-1): f[i] = 0.5*(u[i+1] + u[i]) - 0.5*cfl*(u[i+1] - u[i]) f[n-1] = f[n-2] def BeamWarming_advection(f, u, cfl, n): for i in range(1,n): f[i] = u[i] + 0.5*(1.0 - cfl)*(u[i] - u[i-1]) f[0] = f[1] #------------------------------------------------------------------- # Initial data routines #------------------------------------------------------------------- def id_step(x, args): ul = args[0] ur = args[1] return np.where (x < 0, ul, ur) def exact_step(x, x0, args): ul = args[0] ur = args[1] return np.where(x < x0, ul, ur) #------------------------------------------------------------------- # # id_gaussian: Gaussian initial data # Here args = [A, x0, sigma] # #------------------------------------------------------------------- def id_gaussian(x, args): A = args[0] x0 = args[1] sigma = args[2] return A*np.exp(-(x-x0)**2/sigma**2) def id_box(x, args): ul = args[0] ur = args[1] return np.where (abs(x) < 0.5, ul, ur) def set_bc(f, bctype): N = np.size(f) if (bctype == 1): print "applying outflow bc's. N = " + str(N) f[0] = f[2] f[1] = f[2] f[N-2] = f[N-3] f[N-1] = f[N-3] elif (bctype ==2): f[0] = f[N-4] f[1] = f[N-3] f[N-2] = f[2] f[N-1] = f[3] else: print "unknown bctype = " + str(bctype) exit(2)