Mercurial > repos > public > sbplib

comparison +rv/+time/RungekuttaExteriorRV.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 4b42999874c0
children 5359a61cb4d9
comparison
equal deleted inserted replaced
1016:4b42999874c0 1017:2d7c1333bd6c
4 k % Time step 4 k % Time step
5 t % Time point 5 t % Time point
6 v % Solution vector 6 v % Solution vector
7 n % Time level 7 n % Time level
8 coeffs % The coefficents used for the RK time integration 8 coeffs % The coefficents used for the RK time integration
9 RV % Residual Viscosity 9
10 v_prev % Solution vector at previous time levels, used for the RV update 10 % Properties related to the residual viscositys
11 DvDt % Function for computing the time deriative used for the RV update 11 RV % Residual Viscosity operator
12 lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV update. 12 viscosity % Viscosity vector
13 % dictates which accuracy the boot-strapping should start from. 13 v_prev % Solution vector at previous time levels, used for the RV evaluation
14 upperBdfOrder % Orders of the approximation of the time deriative, used for the RV update. 14 DvDt % Function for computing the time deriative used for the RV evaluation
15 % Dictates the order of accuracy used once the boot-strapping is complete. 15 lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation.
16 % dictates which accuracy the boot-strapping should start from.
17 upperBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation.
18 % Dictates the order of accuracy used once the boot-strapping is complete.
19
20 % Convenience properties. Only for plotting
21 residual
22 dvdt
23 Df
16 end 24 end
17 methods 25 methods
18 26
19 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, rkOrder, bdfOrders) 27 % TODO: Decide on how to compute dvdt.
28 function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, DvDt, rkOrder, bdfOrders)
20 obj.F = F; 29 obj.F = F;
21 obj.k = k; 30 obj.k = k;
22 obj.t = t0; 31 obj.t = t0;
23 obj.v = v0; 32 obj.v = v0;
24 obj.n = 0; 33 obj.n = 0;
32 % the BDF should be done here, or if it should be checked for each 41 % the BDF should be done here, or if it should be checked for each
33 % step taken. 42 % step taken.
34 % If it is moved here, then multiple branching stages can be removed in step() 43 % If it is moved here, then multiple branching stages can be removed in step()
35 % but this will effectively result in a plotted simulation starting from n = upperBdfOrder. 44 % but this will effectively result in a plotted simulation starting from n = upperBdfOrder.
36 % In addition, the properties lowerBdfOrder and upperBdfOrder can be removed. 45 % In addition, the properties lowerBdfOrder and upperBdfOrder can be removed.
37 obj.lowerBdfOrder = bdfOrders.lowerBdfOrder; 46 % obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
38 obj.upperBdfOrder = bdfOrders.upperBdfOrder; 47 % obj.upperBdfOrder = bdfOrders.upperBdfOrder;
39 assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6)); 48 % assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
40 obj.v_prev = []; 49 % obj.v_prev = [];
41 obj.DvDt = rv.time.BdfDerivative(); 50 % obj.DvDt = rv.time.BdfDerivative();
51 % obj.viscosity = zeros(size(v0));
52 % obj.residual = zeros(size(v0));
53 % obj.dvdt = zeros(size(v0));
54 % obj.Df = zeros(size(v0));
55
56 % Using the ODE:
57 obj.DvDt = DvDt;
58 obj.dvdt = obj.DvDt(obj.v);
59 [obj.viscosity, obj.Df] = RV.evaluate(obj.v,obj.dvdt);
60 obj.residual = obj.dvdt + obj.Df;
42 end 61 end
43 62
44 function [v, t] = getV(obj) 63 function [v, t] = getV(obj)
45 v = obj.v; 64 v = obj.v;
46 t = obj.t; 65 t = obj.t;
47 end 66 end
48 67
49 function state = getState(obj) 68 function state = getState(obj)
50 [residual, u_t, grad_f] = obj.RV.getResidual(); 69 state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t);
51 state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
52 end 70 end
53 71
54 function obj = step(obj) 72 function obj = step(obj)
55 % Store current time level 73 % Store current time level and update v_prev
56 numStoredStages = size(obj.v_prev,2); 74 % numStoredStages = size(obj.v_prev,2);
57 if (numStoredStages < obj.upperBdfOrder) 75 % if (numStoredStages < obj.upperBdfOrder)
58 obj.v_prev = [obj.v, obj.v_prev]; 76 % obj.v_prev = [obj.v, obj.v_prev];
59 numStoredStages = numStoredStages+1; 77 % numStoredStages = numStoredStages+1;
60 else 78 % else
61 obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); 79 % obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
62 obj.v_prev(:,1) = obj.v; 80 % obj.v_prev(:,1) = obj.v;
63 end 81 % end
82
83 obj.dvdt = obj.DvDt(obj.v);
84 [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
85 obj.residual = obj.dvdt + obj.Df;
86
64 % Fix the viscosity of the RHS function F 87 % Fix the viscosity of the RHS function F
65 F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV? 88 F_visc = @(v,t) obj.F(v,t,obj.viscosity);
66 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs); 89 obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
67 obj.t = obj.t + obj.k; 90 obj.t = obj.t + obj.k;
68 obj.n = obj.n + 1; 91 obj.n = obj.n + 1;
69 % Calculate dvdt and update RV for the new time level 92
70 if ((numStoredStages >= obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder)) 93 % %Calculate dvdt and evaluate RV for the new time level
71 dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k); 94 % if ((numStoredStages >= obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder))
72 obj.RV.update(obj.v,dvdt); 95 % obj.dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
73 end 96 % [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt);
97 % obj.residual = obj.dvdt + obj.Df;
98 % end
74 end 99 end
75 end 100 end
76 end 101 end