- Thank you received: 1
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
How to implement Crank-Nicolson method: update f^n, f^(n+1)
3 years 9 months ago - 3 years 9 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: 3 years 9 months ago by dong.
3 years 9 months ago - 3 years 9 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: 3 years 9 months ago by hvwahl.
The following user(s) said Thank You: dong
Time to create page: 0.110 seconds
