- Thank you received: 1
How to implement Crank-Nicolson method: update f^n, f^(n+1)
4 years 4 months ago - 4 years 4 months ago #2928
by dong
How to implement Crank-Nicolson method: update f^n, f^(n+1) was created by dong
I tried to implement the Crank-Nicolson method for a transport equation.
The Crank-Nicolson method
[tex]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)[/tex]
or in an incremental form:
[tex]M^*(C^{k+1}-C^k)=-\Delta tAC^k+0.5\Delta t(f^{k+1}+f^k)[/tex]
where [tex]M^*=M+0.5\Delta t A[/tex]
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.
Could you please tell me how to handle this in NGSolve? Thank you so much.
The Crank-Nicolson method
[tex]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)[/tex]
or in an incremental form:
[tex]M^*(C^{k+1}-C^k)=-\Delta tAC^k+0.5\Delta t(f^{k+1}+f^k)[/tex]
where [tex]M^*=M+0.5\Delta t A[/tex]
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.
Code:
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:
Last edit: 4 years 4 months ago by dong.
4 years 4 months ago - 4 years 4 months ago #2942
by hvwahl
Replied by hvwahl on topic How to implement Crank-Nicolson method: update f^n, f^(n+1)
Hi dong,
if you need the right hand side from the last time-step, you can simply store this in a vector
and update this at the end of every time-step:
Best wishes,
Henry
if you need the right hand side from the last time-step, you can simply store this in a vector
Code:
f_last = gfu.vec.CreateVector()
Code:
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
Last edit: 4 years 4 months ago by hvwahl.
The following user(s) said Thank You: dong
Time to create page: 0.112 seconds