Test RK4

src/bvpy/scheme/test.py

 from fenics import *
+import ufl
 import time
 
-T = 2.0            # final time
-num_steps = 50     # number of time steps
-dt = T / num_steps # time step size
+# T = 0.01           # final time
+    # number of time steps
+# dt = T / num_steps # time step size
 
 # Create mesh and define function space
-nx = ny = 30
+nx = ny = 100
 mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
 V = FunctionSpace(mesh, 'P', 1)
 
@@ -26,30 +27,100 @@ u = TrialFunction(V)
 v = TestFunction(V)
 f = Constant(0)
 
-def func(t, state):
-    F = dot(grad(state), grad(v))*dx - f*v*dx
-    return F
+# Create VTK file for saving solutio
+# Time-stepping
 
 
-# Create VTK file for saving solution
-vtkfile = File('heat_gaussian/solution.pvd')
+sol = Function(V)
+sol.assign(u_n)
+k1 = Function(V)
+k2 = Function(V)
+k3 = Function(V)
+k4 = Function(V)
 
-# Time-stepping
+F = (dot(grad(sol), grad(v)))*dx
+
+#Midpoint method
+
+# sol.assign(u_n)
+# scheme = BackwardEuler(F, sol, t=fe.Constant(0), bcs=[bc])
+# solver = RKSolver(scheme)
+# solver.step_interval(0,2,0.04)
+# with XDMFFile('solution.xdmf') as f:
+#     f.write(scheme.solution())
 
 t = 0
-form = (u-u_n)/dt*dx - func(t+0.5*dt, u_n + dt*func(t, u_n))
-a, L = lhs(F), rhs(F)
-u = Function(V)
+vtkfile = XDMFFile('solution.xdmf')
+# vtkfile.write_checkpoint(sol, "heat", t)
+# vtkfile << (sol, t)
+dt = 0.04
+num_steps = 50
 for n in range(num_steps):
-
+    print(t, n)
     # Update current time
-    t += dt
+    vtkfile.write(sol, t)
+
+    evalargs = {u: u_n}
+    F1 = ufl.replace(F, evalargs)
+    # print(F1)
+    assemble(F1, tensor=k1.vector())
+    bc.apply(k1.vector())
+
+    evalargs = {u: u_n + dt*0.5*k1}
+    F2 = ufl.replace(F, evalargs)
+    # print(F2)
+    assemble(F2, tensor=k2.vector())
+    bc.apply(k2.vector())
 
-    # Compute solution
-    solve(a == L, u, bc)
+    evalargs = {u: u_n + dt*0.5*k2}
+    F3 = ufl.replace(F, evalargs)
+    # print(F3)
+    assemble(F3, tensor=k3.vector())
+    bc.apply(k3.vector())
+
+    evalargs = {u: u_n + dt*k3}
+    F4 = ufl.replace(F, evalargs)
+    # print(F4)
+    assemble(F4, tensor=k4.vector())
+    bc.apply(k4.vector())
+
+    form = (u-u_n)/dt*v*dx - ((1/6)*k1+(1/3)*k2+(1/3)*k3+(1/6)*k4)*v*dx
+    # form = (u-u_n)/dt*v*dx - F1
+    a = lhs(form)
+    L = rhs(form)
+    solve(a==L, sol, bc)
 
-    # Save to file and plot solution
-    vtkfile << (u, t)
 
     # Update previous solution
-    u_n.assign(u)
+    u_n.assign(sol)
+    t += dt
+
+vtkfile.write(sol,t)
+
+vtkfile.close()
+
+
+
+# for n in range(num_steps):
+#     print(t, n)
+#     # Update current time
+#
+#     evalargs = {u: u_n}
+#     F1 = ufl.replace(F, evalargs)
+#     assemble(F1, tensor=k1.vector())
+#
+#     evalargs = {u: u_n + dt*0.5*k1}
+#     F2 = ufl.replace(F, evalargs)
+#     assemble(F2, tensor=k2.vector())
+#
+#     form = (u-u_n)/dt*v*dx - k2*v*dx
+#     a = lhs(form)
+#     L = rhs(form)
+#
+#     solve(a==L, sol, bc)
+#
+#     vtkfile << (sol, t)
+#
+#     # Update previous solution
+#     u_n.assign(sol)
+#     t += dt