diff --git a/examples/boussinesq_2d_imex/ProblemClass.py b/examples/boussinesq_2d_imex/ProblemClass.py index 0a9abfa0a1..8ae287234d 100644 --- a/examples/boussinesq_2d_imex/ProblemClass.py +++ b/examples/boussinesq_2d_imex/ProblemClass.py @@ -62,6 +62,7 @@ def __init__(self, cparams, dtype_u, dtype_f): assert 'Nfreq' in cparams assert 'x_bounds' in cparams assert 'z_bounds' in cparams + assert 'order_upw' in cparams assert 'order' in cparams assert 'gmres_maxiter' in cparams assert 'gmres_restart' in cparams @@ -82,7 +83,7 @@ def __init__(self, cparams, dtype_u, dtype_f): self.xx, self.zz, self.h = get2DMesh(self.N, self.x_bounds, self.z_bounds, self.bc_hor[0], self.bc_ver[0]) self.Id, self.M = getBoussinesq2DMatrix(self.N, self.h, self.bc_hor, self.bc_ver, self.c_s, self.Nfreq, self.order) - self.D_upwind = getBoussinesq2DUpwindMatrix( self.N, self.h[0], self.u_adv , self.order) + self.D_upwind = getBoussinesq2DUpwindMatrix( self.N, self.h[0], self.u_adv , self.order_upw) self.logger = logging() diff --git a/examples/boussinesq_2d_imex/build2DFDMatrix.py b/examples/boussinesq_2d_imex/build2DFDMatrix.py index 870a49b9b6..ab8b18ff9f 100644 --- a/examples/boussinesq_2d_imex/build2DFDMatrix.py +++ b/examples/boussinesq_2d_imex/build2DFDMatrix.py @@ -4,11 +4,17 @@ import scipy.sparse as sp from buildFDMatrix import getMatrix, getUpwindMatrix, getBCLeft, getBCRight -def get2DUpwindMatrix(N, dx): +# +# +# +def get2DUpwindMatrix(N, dx, order): - Dx = getUpwindMatrix( N[0], dx) + Dx = getUpwindMatrix( N[0], dx, order) return sp.kron( Dx, sp.eye(N[1]), format="csr" ) +# +# +# def get2DMesh(N, x_b, z_b, bc_hor, bc_ver): assert np.size(N)==2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' assert np.size(x_b)==2, 'x_b needs to be an array with two entries: x_b[0] = left boundary, x_b[1] = right boundary' @@ -39,12 +45,16 @@ def get2DMesh(N, x_b, z_b, bc_hor, bc_ver): xx, zz = np.meshgrid(x,z,indexing="ij") return xx, zz, h -def get2DMatrix(N, h, bc_hor, bc_ver): +# +# +# +def get2DMatrix(N, h, bc_hor, bc_ver, order): + assert np.size(N)==2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz' assert np.size(h)==2, 'h needs to be an array with two entries: h[0]=dx and h[1]=dz' - Ax = getMatrix( N[0], h[0], bc_hor[0], bc_hor[1]) - Az = getMatrix( N[1], h[1], bc_ver[0], bc_ver[1]) + Ax = getMatrix( N[0], h[0], bc_hor[0], bc_hor[1], order) + Az = getMatrix( N[1], h[1], bc_ver[0], bc_ver[1], order) Dx = sp.kron( Ax, sp.eye(N[1]), format="csr") Dz = sp.kron( sp.eye(N[0]), Az, format="csr") @@ -83,4 +93,4 @@ def getBCVertical( value, N, dz, bc_ver): bu = getBCRight( value[1], N[1], dz, bc_ver[1]) bu = np.kron( np.ones(N[0]), bu) - return bd, bu \ No newline at end of file + return bd, bu diff --git a/examples/boussinesq_2d_imex/buildBoussinesq2DMatrix.py b/examples/boussinesq_2d_imex/buildBoussinesq2DMatrix.py index 00768b3991..ce4359c77e 100644 --- a/examples/boussinesq_2d_imex/buildBoussinesq2DMatrix.py +++ b/examples/boussinesq_2d_imex/buildBoussinesq2DMatrix.py @@ -4,7 +4,7 @@ def getBoussinesq2DUpwindMatrix(N, dx, u_adv, order): - Dx = get2DUpwindMatrix(N, dx) + Dx = get2DUpwindMatrix(N, dx, order) # Note: In the equations it is u_t + u_adv* D_x u = ... so in order to comply with the form u_t = M u, # add a minus sign in front of u_adv @@ -19,10 +19,12 @@ def getBoussinesq2DUpwindMatrix(N, dx, u_adv, order): return sp.csc_matrix(M) def getBoussinesq2DMatrix(N, h, bc_hor, bc_ver, c_s, Nfreq, order): - Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0]) - Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1]) - Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2]) - Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3]) + + Dx_u, Dz_u = get2DMatrix(N, h, bc_hor[0], bc_ver[0], order) + Dx_w, Dz_w = get2DMatrix(N, h, bc_hor[1], bc_ver[1], order) + Dx_b, Dz_b = get2DMatrix(N, h, bc_hor[2], bc_ver[2], order) + Dx_p, Dz_p = get2DMatrix(N, h, bc_hor[3], bc_ver[3], order) + Id_N = sp.eye(N[0]*N[1]) Zero = np.zeros((N[0]*N[1],N[0]*N[1])) @@ -53,4 +55,4 @@ def getBoussinesqBCHorizontal(value, N, dx, bc_hor): return b_left, b_right def getBoussinesqBCVertical(): - return 0.0 \ No newline at end of file + return 0.0 diff --git a/examples/boussinesq_2d_imex/buildFDMatrix.py b/examples/boussinesq_2d_imex/buildFDMatrix.py index f42f0d20fa..bbfe7457ac 100644 --- a/examples/boussinesq_2d_imex/buildFDMatrix.py +++ b/examples/boussinesq_2d_imex/buildFDMatrix.py @@ -1,29 +1,38 @@ import numpy as np import scipy.sparse as sp import scipy.linalg as la +import sys + # Only for periodic BC because we have advection only in x direction -def getUpwindMatrix(N, dx): - - #stencil = [-1.0, 1.0] - #zero_pos = 2 - #coeff = 1.0 - - #stencil = [1.0, -4.0, 3.0] - #coeff = 1.0/2.0 - #zero_pos = 3 - - #stencil = [1.0, -6.0, 3.0, 2.0] - #coeff = 1.0/6.0 - #zero_pos = 3 +def getUpwindMatrix(N, dx, order): - #stencil = [-5.0, 30.0, -90.0, 50.0, 15.0] - #coeff = 1.0/60.0 - #zero_pos = 4 + if order==1: + stencil = [-1.0, 1.0] + zero_pos = 2 + coeff = 1.0 - stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0] - coeff = 1.0/60.0 - zero_pos = 5 + elif order==2: + stencil = [1.0, -4.0, 3.0] + coeff = 1.0/2.0 + zero_pos = 3 + + elif order==3: + stencil = [1.0, -6.0, 3.0, 2.0] + coeff = 1.0/6.0 + zero_pos = 3 + + elif order==4: + stencil = [-5.0, 30.0, -90.0, 50.0, 15.0] + coeff = 1.0/60.0 + zero_pos = 4 + elif order==5: + stencil = [3.0, -20.0, 60.0, -120.0, 65.0, 12.0] + coeff = 1.0/60.0 + zero_pos = 5 + else: + sys.exit('Order '+str(order)+' not implemented') + first_col = np.zeros(N) # Because we need to specific first column (not row) in circulant, flip stencil array @@ -34,78 +43,145 @@ def getUpwindMatrix(N, dx): return sp.csc_matrix( coeff*(1.0/dx)*la.circulant(first_col) ) -def getMatrix(N, dx, bc_left, bc_right): - stencil = [1.0, -8.0, 0.0, 8.0, -1.0] - range = [ -2, -1, 0, 1, 2] - A = sp.diags(stencil, range, shape=(N,N)) - A = sp.lil_matrix(A) +def getMatrix(N, dx, bc_left, bc_right, order): assert bc_left in ['periodic','neumann','dirichlet'], "Unknown type of BC" + + if order==2: + stencil = [-1.0, 0.0, 1.0] + range = [-1, 0, 1] + coeff = 1.0/2.0 + elif order==4: + stencil = [1.0, -8.0, 0.0, 8.0, -1.0] + range = [ -2, -1, 0, 1, 2] + coeff = 1.0/12.0 + else: + sys.exit('Order '+str(order)+' not implemented') + + A = sp.diags(stencil, range, shape=(N,N)) + A = sp.lil_matrix(A) + + # + # Periodic boundary conditions + # if bc_left in ['periodic']: assert bc_right in ['periodic'], "Periodic BC can only be selected for both sides simultaneously" if bc_left in ['periodic']: - A[0,N-2] = stencil[0] - A[0,N-1] = stencil[1] - A[1,N-1] = stencil[0] + if order==2: + A[0,N-1] = stencil[0] - if bc_right in ['periodic']: - A[N-2,0] = stencil[4] - A[N-1,0] = stencil[3] - A[N-1,1] = stencil[4] + elif order==4: + A[0,N-2] = stencil[0] + A[0,N-1] = stencil[1] + A[1,N-1] = stencil[0] + if bc_right in ['periodic']: + if order==2: + A[N-1,0] = stencil[2] + elif order==4: + A[N-2,0] = stencil[4] + A[N-1,0] = stencil[3] + A[N-1,1] = stencil[4] + + # + # Neumann boundary conditions + # if bc_left in ['neumann']: A[0,:] = np.zeros(N) - A[0,0] = -8.0 - A[0,1] = 8.0 - A[1,0] = -8.0 + 4.0/3.0 - A[1,1] = -1.0/3.0 + if order==2: + A[0,0] = -4.0/3.0 + A[0,1] = 4.0/3.0 + elif order==4: + A[0,0] = -8.0 + A[0,1] = 8.0 + A[1,0] = -8.0 + 4.0/3.0 + A[1,1] = -1.0/3.0 if bc_right in ['neumann']: A[N-1,:] = np.zeros(N) - A[N-2,N-1] = 8.0 - 4.0/3.0 - A[N-2,N-2] = 1.0/3.0 - A[N-1,N-1] = 8.0 - A[N-1,N-2] = -8.0 - + if order==2: + A[N-1,N-2] = -4.0/3.0 + A[N-1,N-1] = 4.0/3.0 + elif order==4: + A[N-2,N-1] = 8.0 - 4.0/3.0 + A[N-2,N-2] = 1.0/3.0 + A[N-1,N-1] = 8.0 + A[N-1,N-2] = -8.0 + + # + # Dirichlet boundary conditions + # if bc_left in ['dirichlet']: - A[0,:] = np.zeros(N) - A[0,1] = 6.0 + # For order==2, nothing to do here + if order==4: + A[0,:] = np.zeros(N) + A[0,1] = 6.0 if bc_right in ['dirichlet']: - A[N-1,:] = np.zeros(N) - A[N-1,N-2] = -6.0 - - A = 1.0/(12.0*dx)*A + # For order==2, nothing to do here + if order==4: + A[N-1,:] = np.zeros(N) + A[N-1,N-2] = -6.0 + + A = coeff*(1.0/dx)*A return sp.csc_matrix(A) -def getBCLeft(value, N, dx, type): +# +# +# +def getBCLeft(value, N, dx, type, order): assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC" + if order==2: + coeff = 1.0/2.0 + elif order==4: + coeff = 1.0/12.0 + b = np.zeros(N) if type in ['dirichlet']: - b[0] = -6.0*value - b[1] = 1.0*value + if order==2: + b[0] = -value; + elif order==4: + b[0] = -6.0*value + b[1] = 1.0*value if type in ['neumann']: - b[0] = 4.0*dx*value - b[1] = -(2.0/3.0)*dx*value + if order==2: + b[0] = (2.0/3.0)*dx*value + elif order==4: + b[0] = 4.0*dx*value + b[1] = -(2.0/3.0)*dx*value - return (1.0/(12.0*dx))*b + return coeff*(1.0/dx)*b -def getBCRight(value, N, dx, type): +# +# +# +def getBCRight(value, N, dx, type, order): assert type in ['periodic','neumann','dirichlet'], "Unknown type of BC" + if order==2: + coeff = 1.0/2.0 + elif order==4: + coeff = 1.0/12.0 + b = np.zeros(N) if type in ['dirichlet']: - b[N-2] = -1.0*value - b[N-1] = 6.0*value + if order==2: + b[N-1] = value + elif order==4: + b[N-2] = -1.0*value + b[N-1] = 6.0*value if type in ['neumann']: - b[N-2] = -(2.0/3.0)*dx*value - b[N-1] = 4.0*dx*value + if order==2: + b[N-1] = (2.0/3.0)*dx*value + elif order==4: + b[N-2] = -(2.0/3.0)*dx*value + b[N-1] = 4.0*dx*value - return (1.0/(12.0*dx))*b + return coeff*(1.0/dx)*b diff --git a/examples/boussinesq_2d_imex/playground.py b/examples/boussinesq_2d_imex/playground.py index 377c67f407..f75816e03b 100644 --- a/examples/boussinesq_2d_imex/playground.py +++ b/examples/boussinesq_2d_imex/playground.py @@ -24,11 +24,11 @@ # set global logger (remove this if you do not want the output at all) logger = Log.setup_custom_logger('root') - num_procs = 1 + num_procs = 8 # This comes as read-in for the level class lparams = {} - lparams['restol'] = 1E-14 + lparams['restol'] = 1E-8 swparams = {} swparams['collocation_class'] = collclass.CollGaussLobatto @@ -36,28 +36,29 @@ swparams['do_LU'] = True sparams = {} - sparams['maxiter'] = 4 + sparams['maxiter'] = 12 # setup parameters "in time" t0 = 0 Tend = 3000 Nsteps = 500 - Tend = 30 - Nsteps = 5 + #Tend = 30 + #Nsteps = 5 dt = Tend/float(Nsteps) # This comes as read-in for the problem class pparams = {} - pparams['nvars'] = [(4,450,30)] + pparams['nvars'] = [(4,450,30), (4,450,30)] pparams['u_adv'] = 0.02 pparams['c_s'] = 0.3 pparams['Nfreq'] = 0.01 pparams['x_bounds'] = [(-150.0, 150.0)] pparams['z_bounds'] = [( 0.0, 10.0)] - pparams['order'] = [0] # [fine_level, coarse_level] - pparams['gmres_maxiter'] = [50] - pparams['gmres_restart'] = 20 - pparams['gmres_tol'] = 1e-6 + pparams['order'] = [4, 2] # [fine_level, coarse_level] + pparams['order_upw'] = [5, 1] + pparams['gmres_maxiter'] = [50, 50] + pparams['gmres_restart'] = [20, 20] + pparams['gmres_tol'] = [1e-8, 1e-8] # This comes as read-in for the transfer operations tparams = {} @@ -73,8 +74,8 @@ description['sweeper_class'] = imex_1st_order description['level_params'] = lparams description['hook_class'] = plot_solution - #description['transfer_class'] = mesh_to_mesh_2d - #description['transfer_params'] = tparams + description['transfer_class'] = mesh_to_mesh_2d + description['transfer_params'] = tparams # quickly generate block of steps MS = mp.generate_steps(num_procs,sparams,description) @@ -105,4 +106,4 @@ #extract_stats = grep_stats(stats,iter=-1,type='residual') #sortedlist_stats = sort_stats(extract_stats,sortby='step') - #print(extract_stats,sortedlist_stats) \ No newline at end of file + #print(extract_stats,sortedlist_stats)