Mercurial > repos > public > sbplib

comparison +rv/+time/rungekuttaRV.m @ 1017:2d7c1333bd6c feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using the ODE to approximate the time derivative in the residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Dec 2018 16:29:21 +0100
parents 1e437c9e5132
children 010bb2677230
comparison
equal deleted inserted replaced
1016:4b42999874c0 1017:2d7c1333bd6c
1 % Takes one time step of size dt using the rungekutta method 1 % Takes one time step of size dt using the rungekutta method
2 % starting from v and where the function F(v,t,RV) gives the 2 % starting from v and where the function F(v,t,RV) gives the
3 % time derivatives. coeffs is a struct holding the RK coefficients 3 % time derivatives. coeffs is a struct holding the RK coefficients
4 % for the specific method. RV is the residual viscosity which is updated 4 % for the specific method. RV is the residual viscosity which is updated
5 % in between the stages and after the updated solution is computed. 5 % in between the stages and after the updated solution is computed.
6 function v = rungekuttaRV(v, t , dt, F, RV, coeffs) 6 function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs)
7 % Move one stage outside to avoid branching for updating the 7 % Move one stage outside to avoid branching for updating the
8 % residual inside the loop. 8 % residual inside the loop.
9 k = zeros(length(v), coeffs.s); 9 k = zeros(length(v), coeffs.s);
10 k(:,1) = F(v,t,RV.getViscosity()); 10 k(:,1) = F(v,t,RV.evaluate(v,DvDt(v)));
11 11
12 % Compute the intermediate stages k 12 % Compute the intermediate stages k
13 for i = 2:coeffs.s 13 for i = 2:coeffs.s
14 u = v; 14 u = v;
15 for j = 1:i-1 15 for j = 1:i-1
16 u = u + dt*coeffs.a(i,j)*k(:,j); 16 u = u + dt*coeffs.a(i,j)*k(:,j);
17 end 17 end
18 RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization 18 %RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization
19 k(:,i) = F(u,t+coeffs.c(i)*dt, RV.getViscosity()); 19 k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluate(u,DvDt(u)));
20 end 20 end
21 21
22 % Compute the updated solution as a linear combination 22 % Compute the updated solution as a linear combination
23 % of the intermediate stages. 23 % of the intermediate stages.
24 u = v; 24 u = v;
25 for i = 1:coeffs.s 25 for i = 1:coeffs.s
26 u = u + dt*coeffs.b(i)*k(:,i); 26 u = u + dt*coeffs.b(i)*k(:,i);
27 end 27 end
28 RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization 28 %RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization
29 v = u; 29 v = u;
30 end 30 end