Mercurial > repos > public > sbplib

changeset 835:008496ca38f3 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Compute the residual in between each runge-kutta stage. Note: It is not clear whether the correct residual is used when computing the stages. Must investigate further.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 13 Sep 2018 18:14:54 +0200
parents f1f0bf087e1c
children 619561e9ec0e
files +time/+rk4/rungekutta_4RV.m +time/Rungekutta4RV.m
diffstat 2 files changed, 22 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk4/rungekutta_4RV.m	Thu Sep 13 18:14:54 2018 +0200
@@ -0,0 +1,18 @@
+% Takes one time step of size k using the rungekutta method
+% starting from v_0 and where the function F(v,t) gives the
+% time derivatives.
+function v = rungekutta_4RV(v, t , k, F, RV)
+    k1 = F(v, t, RV.getViscosity());
+	
+    RV.update(v+0.5*k*k1, v, 0.5*k);
+	k2 = F(v+0.5*k*k1, t+0.5*k, RV.getViscosity());
+	
+	RV.update(v+0.5*k*k2, v, 0.5*k);
+    k3 = F(v+0.5*k*k2, t+0.5*k, RV.getViscosity());
+    
+    RV.update(v+k*k3, v, k);
+    k4 = F(v+k*k3,t+k, RV.getViscosity());
+    
+    RV.update(v + (1/6)*(k1+2*(k2+k3)+k4)*k, v, k);
+    v = v + (1/6)*(k1+2*(k2+k3)+k4)*k;
+end
\ No newline at end of file
--- a/+time/Rungekutta4RV.m	Thu Sep 13 14:21:47 2018 +0200
+++ b/+time/Rungekutta4RV.m	Thu Sep 13 18:14:54 2018 +0200
@@ -28,17 +28,13 @@
             t = obj.t;
         end
 
-        function [residual, viscosity, t] = getRV(obj)
-            residual = obj.RV.getResidual();
-            viscosity = obj.RV.getViscosity();
-            t = obj.t;
+        function state = getState(obj)
+            [residual, u_t, f_x] = obj.RV.getResidual();
+            state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'f_x', f_x, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
         end
 
         function obj = step(obj)
-            v_prev = obj.v;
-            F = @(v,t) obj.F(v, t, obj.RV.getViscosity());
-            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, F);
-            obj.RV.update(obj.v, v_prev, obj.k);
+            obj.v = time.rk4.rungekutta_4RV(obj.v, obj.t, obj.k, obj.F, obj.RV);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end