diff --git a/examples/boussinesq_2d_imex/plotgmrescounter.py b/examples/boussinesq_2d_imex/plotgmrescounter.py index c568adfbde..5ea2a09921 100644 --- a/examples/boussinesq_2d_imex/plotgmrescounter.py +++ b/examples/boussinesq_2d_imex/plotgmrescounter.py @@ -14,7 +14,9 @@ udirk = np.load('dirk.npy') uimex = np.load('rkimex.npy') uref = np.load('uref.npy') + usplit = np.load('split.npy') + print("Estimated discretisation error split explicit: %5.3e" % ( np.linalg.norm(usplit.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )) print("Estimated discretisation error of DIRK: %5.3e" % ( np.linalg.norm(udirk.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )) print("Estimated discretisation error of RK-IMEX: %5.3e" % ( np.linalg.norm(uimex.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )) print("Estimated discretisation error of SDC: %5.3e" % ( np.linalg.norm(uend.flatten() - uref.flatten(), np.inf)/np.linalg.norm(uref.flatten(),np.inf) )) diff --git a/examples/boussinesq_2d_imex/rungmrescounter.py b/examples/boussinesq_2d_imex/rungmrescounter.py index e7a7cd295b..1b6da42c4c 100644 --- a/examples/boussinesq_2d_imex/rungmrescounter.py +++ b/examples/boussinesq_2d_imex/rungmrescounter.py @@ -21,7 +21,7 @@ from unflatten import unflatten -from standard_integrators import dirk, bdf2, trapezoidal, rk_imex +from standard_integrators import dirk, bdf2, trapezoidal, rk_imex, SplitExplicit if __name__ == "__main__": @@ -46,8 +46,8 @@ # setup parameters "in time" t0 = 0 - Tend = 3000 - Nsteps = 500 + Tend = 3000 + Nsteps = 500 dt = Tend/float(Nsteps) # This comes as read-in for the problem class @@ -96,19 +96,33 @@ print("CFL number of acoustics (horizontal): %4.2f" % cfl_acoustic_hor) print("CFL number of acoustics (vertical): %4.2f" % cfl_acoustic_ver) - dirkp = dirk(P, np.min([4,dirk_order])) + method_split = 'MIS4_4' +# method_split = 'RK3' + splitp = SplitExplicit(P, method_split, pparams) u0 = uinit.values.flatten() + usplit = np.copy(u0) +# for i in range(0,Nsteps): + print(np.linalg.norm(usplit)) + for i in range(0,2*Nsteps): + usplit = splitp.timestep(usplit, dt/2) + print(np.linalg.norm(usplit)) + + dirkp = dirk(P, np.min([4,dirk_order])) udirk = np.copy(u0) print("Running DIRK ....") + print(np.linalg.norm(udirk)) for i in range(0,Nsteps): udirk = dirkp.timestep(udirk, dt) + print np.linalg.norm(udirk) rkimex = rk_imex(P, dirk_order) uimex = np.copy(u0) dt_imex = dt print("Running RK-IMEX ....") for i in range(0,Nsteps): +# print("Running RK-IMEWX Step: %4.2f" % dt_imex) uimex = rkimex.timestep(uimex, dt_imex) + print(np.linalg.norm(uimex)) # call main function to get things done... print("Running SDC...") @@ -123,6 +137,7 @@ for i in range(0,10*Nsteps): uref = rkimexref.timestep(uref, dt_ref) + usplit = unflatten(usplit, 4, P.N[0], P.N[1]) udirk = unflatten(udirk, 4, P.N[0], P.N[1]) uimex = unflatten(uimex, 4, P.N[0], P.N[1]) uref = unflatten(uref, 4, P.N[0], P.N[1]) @@ -131,8 +146,17 @@ np.save('sdc', uend.values) np.save('dirk', udirk) np.save('rkimex', uimex) + np.save('split', usplit) np.save('uref', uref) + + print "diff split ",np.linalg.norm(uref-usplit) + print "diff dirk ",np.linalg.norm(uref-udirk) + print "diff rkimex ",np.linalg.norm(uref-uimex) + print "diff sdc ",np.linalg.norm(uref-uend.values) + print(" #### Logging report for Split #### " ) + print("Total number of matrix multiplcations: %5i" % splitp.logger.nsmall) + print(" #### Logging report for DIRK-%1i #### " % dirkp.order) print("Number of calls to implicit solver: %5i" % dirkp.logger.solver_calls) print("Total number of GMRES iterations: %5i" % dirkp.logger.iterations) diff --git a/examples/boussinesq_2d_imex/standard_integrators.py b/examples/boussinesq_2d_imex/standard_integrators.py index 339a1fe65b..16a0796bb5 100644 --- a/examples/boussinesq_2d_imex/standard_integrators.py +++ b/examples/boussinesq_2d_imex/standard_integrators.py @@ -270,8 +270,149 @@ def f_solve(self, b, alpha, u0): return sol # -# A diagonally implicit Runge-Kutta method of order 2, 3 or 4 +# Split-Explicit method # + +class SplitExplicit: + + def __init__(self, problem, method, pparams): + + assert isinstance(problem, boussinesq_2d_imex), "problem is wrong type of object" + self.Ndof = np.shape(problem.M)[0] + self.method = method + self.logger = logging() + self.problem = problem + self.pparams = pparams + self.NdofTher = 2*problem.N[0]*problem.N[1] + self.NdofMom = 2*problem.N[0]*problem.N[1] + + print "dx ",problem.h[0] + print "dz ",problem.h[1] + + assert self.method in ["MIS4_4","RK3"], 'Method must be MIS4_4' + + if (self.method=='RK3'): + self.nstages=3 + self.aRunge=np.zeros((4,4)) + self.aRunge[0,0]= 1./3. + self.aRunge[1,1]= 1./2. + self.aRunge[2,2]= 1. + self.dRunge=np.zeros((4,4)) + self.gRunge=np.zeros((4,4)) + if (self.method=='MIS4_4'): + self.nstages=4 + self.aRunge=np.zeros((4,4)) + self.aRunge[0,0]= 0.38758444641450318 + self.aRunge[1,0]= -2.5318448354142823e-002 + self.aRunge[1,1]= 0.38668943087310403 + self.aRunge[2,0]= 0.20899983523553325 + self.aRunge[2,1]= -0.45856648476371231 + self.aRunge[2,2]= 0.43423187573425748 + self.aRunge[3,0]= -0.10048822195663100 + self.aRunge[3,1]= -0.46186171956333327 + self.aRunge[3,2]= 0.83045062122462809 + self.aRunge[3,3]= 0.27014914900250392 + self.dRunge=np.zeros((4,4)) + self.dRunge[1,1]= 0.52349249922385610 + self.dRunge[2,1]= 1.1683374366893629 + self.dRunge[2,2]= -0.75762080241712637 + self.dRunge[3,1]= -3.6477233846797109e-002 + self.dRunge[3,2]= 0.56936148730740477 + self.dRunge[3,3]= 0.47746263002599681 + self.gRunge=np.zeros((4,4)) + self.gRunge[1,1]= 0.13145089796226542 + self.gRunge[2,1]= -0.36855857648747881 + self.gRunge[2,2]= 0.33159232636600550 + self.gRunge[3,1]= -6.5767130537473045E-002 + self.gRunge[3,2]= 4.0591093109036858E-002 + self.gRunge[3,3]= 6.4902111640806712E-002 + self.dtRunge=np.zeros(self.nstages); + for i in range(0,self.nstages): + self.dtRunge[i]=0 + temp=1. + for j in range(0,i+1): + self.dtRunge[i]=self.dtRunge[i]+self.aRunge[i,j] + temp=temp-self.dRunge[i,j] + self.dRunge[i,0]=temp + for j in range(0,i+1): + self.aRunge[i,j]=self.aRunge[i,j]/self.dtRunge[i] + self.gRunge[i,j]=self.gRunge[i,j]/self.dtRunge[i] + + self.U = np.zeros((self.Ndof,self.nstages+1)) + self.F = np.zeros((self.Ndof,self.nstages)) + self.FSlow = np.zeros((self.Ndof)) + self.nsMin = 8 + self.logger.nsmall = 0 + + def NumSmallTimeSteps(self , dx, dz, dt): + + cs = self.pparams['c_s'] + ns = dt / (.9 / np.sqrt(1/(dx*dx)+1/(dz*dz)) / cs) + ns = max(np.int(np.ceil(ns)),self.nsMin) +# print " dx " , dx +# print " dz " , dz +# print " cs " , cs +# print 1/np.sqrt(1/(dx*dx)+1/(dz*dz))/cs +# print " dt ",dt," ns ",ns + + return ns + + + def timestep(self, u0, dt): + + self.U[:,0] = u0 + + self.ns = self.NumSmallTimeSteps(self.problem.h[0], self.problem.h[1], dt) + + for i in range(0,self.nstages): + self.F[:,i] = self.f_slow(self.U[:,i]) + self.FSlow[:] = 0. + for j in range(0,i+1): + self.FSlow += (self.aRunge[i,j]*self.F[:,j] + self.gRunge[i,j]/dt*(self.U[:,j] - u0)) + self.U[:,i+1] = 0 + for j in range(0,i+1): + self.U[:,i+1] += self.dRunge[i,j]*self.U[:,j] + nsLoc = np.int(np.ceil(self.ns*self.dtRunge[i])) + self.logger.nsmall += nsLoc + dtLoc = dt*self.dtRunge[i] + dTau = dtLoc/nsLoc + self.U[:,i+1] = self.VerletLin(self.U[:,i+1], self.FSlow, nsLoc, dTau) + u0 = self.U[:,self.nstages] + return u0 + + + def VerletLin(self, u0, FSlow, ns, dTau): + for i in range(0,ns): + u0[0:self.NdofMom] += dTau * (self.f_fastMom(u0) + FSlow[0:self.NdofMom]) + u0[self.NdofMom:self.Ndof] += dTau * (self.f_fastTher(u0) + FSlow[self.NdofMom:self.Ndof]) + + return u0 + + def RK3Lin(self, u0, FSlow, ns, dTau): + + u = u0 + for i in range(0,ns): + u = u0 + dTau/3.*(self.f_fast(u) + FSlow) + u = u0 + dTau/2.*(self.f_fast(u) + FSlow) + u = u0 + dTau*(self.f_fast(u) + FSlow) + u0 = u + + return u0 + + def f_slow(self, u): + return self.problem.D_upwind.dot(u) + + def f_fast(self, u): + return self.problem.M.dot(u) + + def f_fastMom(self, u): + return self.problem.M[0:self.NdofMom,self.NdofMom:self.Ndof]\ + .dot(u[self.NdofMom:self.Ndof]) + + def f_fastTher(self, u): + return self.problem.M[self.NdofMom:self.Ndof,0:self.NdofMom]\ + .dot(u[0:self.NdofMom]) + class dirk: def __init__(self, problem, order): @@ -283,7 +424,7 @@ def __init__(self, problem, order): self.problem = problem assert self.order in [2,22,3,4], 'Order must be 2,22,3,4' - + if (self.order==2): self.nstages = 1 self.A = np.zeros((1,1))