How to implement Crank-Nicolson method: update f^n, f^(n+1)

1 week 2 days ago - 1 week 2 days ago #2928 by dong
I tried to implement the Crank-Nicolson method for a transport equation.

The Crank-Nicolson method
\[M\frac{C^{k+1}-C^k}{\Delta t}+\frac{1}{2}A[ C^{k+1}+C^k]=\frac{1}{2}(f^{k+1}+f^k)\]
or in an incremental form:
\[M^*(C^{k+1}-C^k)=-\Delta tAC^k+0.5\Delta t(f^{k+1}+f^k)\]
where
\[M^*=M+0.5\Delta t A\]

I don't know how to represent (f^{k+1}+f^k) in the time loop. I searched on the tutorial, but I didn't find an example. Please see the attached file.
while t_intermediate < Tend-0.5*dt:
    #Update time parameter
    t.Set(time+t_intermediate+dt)
    f.Assemble()
    
    #RHS: 1/2*dt*(f^(n+1)+f^(n))-dt*A*C^(k)??
    res.data = 0.5*dt * (f.vec+f.vec) - dt * a.mat * gfu.vec #wrong
    gfu.vec.data += invmstar * res
        
    t_intermediate += dt
    #print("\r", time + t_intermediate,   end="")
    Redraw(blocking=True)

time += t_intermediate

Could you please tell me how to handle this in NGSolve? Thank you so much.
Attachments:

Please Log in or Create an account to join the conversation.

1 week 1 day ago - 1 week 1 day ago #2942 by hvwahl
Hi dong,

if you need the right hand side from the last time-step, you can simply store this in a vector
f_last = gfu.vec.CreateVector()
and update this at the end of every time-step:
f.Assemble()
f_last.data = f.vec 
while t_intermediate < Tend-0.5*dt:
    t.Set(time+t_intermediate+dt)
    f.Assemble()

    res.data = 0.5*dt * (f.vec+f_last) - dt * a.mat * gfu.vec
    gfu.vec.data += invmstar * res
        
    f_last.data = f.vec
    t_intermediate += dt

Best wishes,
Henry
The following user(s) said Thank You: dong

Please Log in or Create an account to join the conversation.

6 days 23 hours ago #2947 by dong
Thank you so much, Henry. I got it.

Please Log in or Create an account to join the conversation.

© 2019 Netgen/NGSolve