1*55a74a43SLisandro Dalcin# Stiff scalar valued ODE problem with an exact solution 2*55a74a43SLisandro Dalcin 3*55a74a43SLisandro Dalcinimport sys, petsc4py 4*55a74a43SLisandro Dalcinpetsc4py.init(sys.argv) 5*55a74a43SLisandro Dalcin 6*55a74a43SLisandro Dalcinfrom petsc4py import PETSc 7*55a74a43SLisandro Dalcinfrom math import sin, cos, exp 8*55a74a43SLisandro Dalcin 9*55a74a43SLisandro Dalcinclass CE(object): 10*55a74a43SLisandro Dalcin n = 1 11*55a74a43SLisandro Dalcin comm = PETSc.COMM_SELF 12*55a74a43SLisandro Dalcin def __init__(self, lambda_=1.0): 13*55a74a43SLisandro Dalcin self.lambda_ = lambda_ 14*55a74a43SLisandro Dalcin def evalSolution(self, t, x): 15*55a74a43SLisandro Dalcin l = self.lambda_ 16*55a74a43SLisandro Dalcin x[0] = l/(l*l+1)*(l*cos(t)+sin(t)) - l*l/(l*l+1)*exp(-l*t) 17*55a74a43SLisandro Dalcin x.assemble() 18*55a74a43SLisandro Dalcin def evalFunction(self, ts, t, x, xdot, f): 19*55a74a43SLisandro Dalcin l = self.lambda_ 20*55a74a43SLisandro Dalcin f[0] = xdot[0] + l*(x[0] - cos(t)) 21*55a74a43SLisandro Dalcin f.assemble() 22*55a74a43SLisandro Dalcin def evalJacobian(self, ts, t, x, xdot, a, A, B): 23*55a74a43SLisandro Dalcin J = B 24*55a74a43SLisandro Dalcin l = self.lambda_ 25*55a74a43SLisandro Dalcin J[0,0] = a + l 26*55a74a43SLisandro Dalcin J.assemble() 27*55a74a43SLisandro Dalcin if A != B: A.assemble() 28*55a74a43SLisandro Dalcin return True # same nonzero pattern 29*55a74a43SLisandro Dalcin 30*55a74a43SLisandro DalcinOptDB = PETSc.Options() 31*55a74a43SLisandro Dalcin 32*55a74a43SLisandro Dalcinlambda_ = OptDB.getScalar('lambda', 10.0) 33*55a74a43SLisandro Dalcinode = CE(lambda_) 34*55a74a43SLisandro Dalcin 35*55a74a43SLisandro DalcinJ = PETSc.Mat().createDense([ode.n,ode.n], comm=ode.comm) 36*55a74a43SLisandro DalcinJ.setUp() 37*55a74a43SLisandro Dalcinx = PETSc.Vec().createSeq(ode.n, comm=ode.comm) 38*55a74a43SLisandro Dalcinf = x.duplicate() 39*55a74a43SLisandro Dalcin 40*55a74a43SLisandro Dalcints = PETSc.TS().create(comm=ode.comm) 41*55a74a43SLisandro Dalcints.setProblemType(ts.ProblemType.NONLINEAR) 42*55a74a43SLisandro Dalcints.setType(ts.Type.GLLE) 43*55a74a43SLisandro Dalcin 44*55a74a43SLisandro Dalcints.setIFunction(ode.evalFunction, f) 45*55a74a43SLisandro Dalcints.setIJacobian(ode.evalJacobian, J) 46*55a74a43SLisandro Dalcin 47*55a74a43SLisandro Dalcints.setTime(0.0) 48*55a74a43SLisandro Dalcints.setTimeStep(0.001) 49*55a74a43SLisandro Dalcints.setMaxTime(10) 50*55a74a43SLisandro Dalcints.setMaxSteps(10000) 51*55a74a43SLisandro Dalcints.setExactFinalTime(PETSc.TS.ExactFinalTime.INTERPOLATE) 52*55a74a43SLisandro Dalcin 53*55a74a43SLisandro Dalcinclass Monitor(object): 54*55a74a43SLisandro Dalcin def __init__(self, ode): 55*55a74a43SLisandro Dalcin self.ode = ode 56*55a74a43SLisandro Dalcin self.x = PETSc.Vec().createSeq(ode.n, comm=ode.comm) 57*55a74a43SLisandro Dalcin def __call__(self, ts, k, t, x): 58*55a74a43SLisandro Dalcin ode.evalSolution(t, self.x) 59*55a74a43SLisandro Dalcin self.x.axpy(-1, x) 60*55a74a43SLisandro Dalcin e = self.x.norm() 61*55a74a43SLisandro Dalcin h = ts.getTimeStep() 62*55a74a43SLisandro Dalcin PETSc.Sys.Print("step %3d t=%8.2e h=%8.2e error=%8.2e" % 63*55a74a43SLisandro Dalcin (k, t, h, e), comm=self.ode.comm) 64*55a74a43SLisandro Dalcin 65*55a74a43SLisandro Dalcints.setMonitor(Monitor(ode)) 66*55a74a43SLisandro Dalcin 67*55a74a43SLisandro Dalcints.setFromOptions() 68*55a74a43SLisandro Dalcinode.evalSolution(0.0, x) 69*55a74a43SLisandro Dalcints.solve(x) 70*55a74a43SLisandro Dalcin 71*55a74a43SLisandro Dalcindel ode, J, x, f, ts 72