diff --git a/src/bvpy/scheme/test.py b/src/bvpy/scheme/test.py index 0504c1c763f329c5cdd9ba6366e400bf4cfe4ab4..a0e926f99a47ad05957d01d812a7202e1d39ef99 100644 --- a/src/bvpy/scheme/test.py +++ b/src/bvpy/scheme/test.py @@ -1,12 +1,13 @@ 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