Mercurial > repos > public > sbplib
changeset 1169:d02e5b8a0b24 feature/rv
Rename RungekuttaRV time steppers. Add RungekuttaRVMultiStage time stepper
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 28 Jun 2019 13:13:17 +0200 |
| parents | af3c4eb0cbbd |
| children | d3bde8a23e08 |
| files | +rv/+time/RungeKuttaRvInstage.m +rv/+time/RungekuttaExteriorRv.m +rv/+time/RungekuttaExteriorRvBdf.m +rv/+time/RungekuttaExteriorRvMg.m +rv/+time/RungekuttaExteriorRvMultiStage.m +rv/+time/RungekuttaInteriorRv.m +rv/+time/RungekuttaRvBdf.m +rv/+time/RungekuttaRvMultiGrid.m |
| diffstat | 8 files changed, 270 insertions(+), 263 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungeKuttaRvInstage.m Fri Jun 28 13:13:17 2019 +0200 @@ -0,0 +1,47 @@ +classdef RungekuttaRvInstage < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + coeffs % The coefficents used for the RK time integration + RV % Residual Viscosity + DvDt % Function for computing the time deriative used for the RV evaluation + end + + methods + function obj = RungekuttaRvInstage(F, k, t0, v0, RV, DvDt, order) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(order); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + obj.RV = RV; + obj.DvDt = DvDt; + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + dvdt = obj.DvDt(obj.v); + [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); + state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); + end + + % Advances the solution vector one time step using the Runge-Kutta method given by + % obj.coeffs, updating the Residual Viscosity in each Runge-Kutta stage + function obj = step(obj) + obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- a/+rv/+time/RungekuttaExteriorRv.m Fri Jun 28 13:12:29 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -classdef RungekuttaExteriorRv < time.Timestepper - properties - F % RHS of the ODE - k % Time step - t % Time point - v % Solution vector - n % Time level - rkScheme % The particular RK scheme used for time integration - RV % Residual Viscosity operator - DvDt % Function for computing the time deriative used for the RV evaluation - end - methods - - function obj = RungekuttaExteriorRv(F, k, t0, v0, RV, DvDt, order) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.n = 0; - - if (order == 4) % Use specialized RK4 scheme - obj.rkScheme = @time.rk.rungekutta_4; - else - % Extract the coefficients for the specified order - % used for the RK updates from the Butcher tableua. - [s,a,b,c] = time.rk.butcherTableau(order); - coeffs = struct('s',s,'a',a,'b',b,'c',c); - obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); - end - - obj.RV = RV; - obj.DvDt = DvDt; - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - dvdt = obj.DvDt(obj.v); - [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); - state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); - end - - % Advances the solution vector one time step using the Runge-Kutta method given by - % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps - function obj = step(obj) - % Fix the viscosity of the RHS function F - viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v)); - m = length(viscosity); - F_visc = @(v,t) obj.F(v,t,spdiags(viscosity,0,m,m)); - obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_visc); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end -end \ No newline at end of file
--- a/+rv/+time/RungekuttaExteriorRvBdf.m Fri Jun 28 13:12:29 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -classdef RungekuttaExteriorRvBdf < time.Timestepper - properties - F % RHS of the ODE - k % Time step - t % Time point - v % Solution vector - n % Time level - rkScheme % The particular RK scheme used for time integration - - - % Properties related to the residual viscositys - RV % Residual Viscosity operator - v_prev % Solution vector at previous time levels, used for the RV evaluation - DvDt % Function for computing the time deriative used for the RV evaluation - lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation. - % dictates which accuracy the boot-strapping should start from. - upperBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation. - % Dictates the order of accuracy used once the boot-strapping is complete. - - - end - methods - function obj = RungekuttaExteriorRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.n = 0; - obj.RV = RV; - obj.lowerBdfOrder = bdfOrders.lowerBdfOrder; - obj.upperBdfOrder = bdfOrders.upperBdfOrder; - assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6)); - obj.v_prev = []; - obj.DvDt = rv.time.BdfDerivative(); - - if (rkOrder == 4) % Use specialized RK4 scheme - obj.rkScheme = @time.rk.rungekutta_4; - else - % Extract the coefficients for the specified order - % used for the RK updates from the Butcher tableua. - [s,a,b,c] = time.rk.butcherTableau(rkOrder); - coeffs = struct('s',s,'a',a,'b',b,'c',c); - obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); - end - - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - if (size(obj.v_prev,2) >= obj.lowerBdfOrder) - dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k); - [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); - else - viscosity = zeros(size(obj.v)); - dvdt = zeros(size(obj.v)); - Df = zeros(size(obj.v)); - firstOrderViscosity = zeros(size(obj.v)); - residualViscosity = zeros(size(obj.v)); - end - state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); - end - - function obj = step(obj) - nStoredStages = size(obj.v_prev,2); - - %Calculate viscosity for the new time level - if (nStoredStages >= obj.lowerBdfOrder) - viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k)); - else - viscosity = zeros(size(obj.v)); - end - - % Store current time level and update v_prev - if (nStoredStages < obj.upperBdfOrder) - obj.v_prev = [obj.v, obj.v_prev]; - else - obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); - obj.v_prev(:,1) = obj.v; - end - - % Fix the viscosity of the RHS function F - m = length(viscosity); - F_visc = @(v,t) obj.F(v, t, spdiags(viscosity,0,m,m)); - obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_visc); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end -end \ No newline at end of file
--- a/+rv/+time/RungekuttaExteriorRvMg.m Fri Jun 28 13:12:29 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,65 +0,0 @@ -classdef RungekuttaExteriorRvMg < time.Timestepper - properties - F % RHS of the ODE - F_unstable % RHS of the unstabalized ODE - k % Time step - t % Time point - v % Solution vector - n % Time level - rkScheme % The particular RK scheme used for time integration - RV % Residual Viscosity operator - DvDt % Function for computing the time deriative used for the RV evaluation - v_unstable - viscosity - end - methods - - function obj = RungekuttaExteriorRvMg(F, F_unstable, k, t0, v0, RV, DvDt, order) - obj.F = F; - obj.F_unstable = F_unstable; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.n = 0; - - if (order == 4) % Use specialized RK4 scheme - obj.rkScheme = @time.rk.rungekutta_4; - else - % Extract the coefficients for the specified order - % used for the RK updates from the Butcher tableua. - [s,a,b,c] = time.rk.butcherTableau(order); - coeffs = struct('s',s,'a',a,'b',b,'c',c); - obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); - end - - obj.RV = RV; - obj.DvDt = DvDt; - obj.v_unstable = 0*v0; - obj.viscosity = 0*v0; - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - dvdt = obj.DvDt(obj.v_unstable); - [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); - state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', obj.viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); - end - - % Advances the solution vector one time step using the Runge-Kutta method given by - % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps - function obj = step(obj) - m = length(obj.viscosity); - obj.v_unstable = obj.rkScheme(obj.v, obj.t, obj.k, obj.F_unstable); - obj.viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v_unstable)); - % Fix the viscosity of the stabilized RHS - F_stable = @(v,t) obj.F(v,t,spdiags(obj.viscosity,0,m,m)); - obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_stable); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaExteriorRvMultiStage.m Fri Jun 28 13:13:17 2019 +0200 @@ -0,0 +1,65 @@ +classdef RungekuttaRvMultiStage < time.Timestepper + properties + F % RHS of the ODE + F_unstable % RHS of the unstabalized ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + rkScheme % The particular RK scheme used for time integration + RV % Residual Viscosity operator + DvDt % Function for computing the time deriative used for the RV evaluation + v_unstable + viscosity + end + methods + + function obj = RungekuttaRvMultiStage(F, F_unstable, k, t0, v0, RV, DvDt, order) + obj.F = F; + obj.F_unstable = F_unstable; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + + if (order == 4) % Use specialized RK4 scheme + obj.rkScheme = @time.rk.rungekutta_4; + else + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(order); + coeffs = struct('s',s,'a',a,'b',b,'c',c); + obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); + end + + obj.RV = RV; + obj.DvDt = DvDt; + obj.v_unstable = 0*v0; + obj.viscosity = 0*v0; + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + dvdt = obj.DvDt(obj.v_unstable); + [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); + state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', obj.viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); + end + + % Advances the solution vector one time step using the Runge-Kutta method given by + % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps + function obj = step(obj) + m = length(obj.viscosity); + obj.v_unstable = obj.rkScheme(obj.v, obj.t, obj.k, obj.F_unstable); + obj.viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v_unstable)); + % Fix the viscosity of the stabilized RHS + F_stable = @(v,t) obj.F(v,t,spdiags(obj.viscosity,0,m,m)); + obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_stable); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- a/+rv/+time/RungekuttaInteriorRv.m Fri Jun 28 13:12:29 2019 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -classdef RungekuttaInteriorRv < time.Timestepper - properties - F % RHS of the ODE - k % Time step - t % Time point - v % Solution vector - n % Time level - coeffs % The coefficents used for the RK time integration - RV % Residual Viscosity - DvDt % Function for computing the time deriative used for the RV evaluation - end - - methods - function obj = RungekuttaInteriorRv(F, k, t0, v0, RV, DvDt, order) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.n = 0; - % Extract the coefficients for the specified order - % used for the RK updates from the Butcher tableua. - [s,a,b,c] = time.rk.butcherTableau(order); - obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); - obj.RV = RV; - obj.DvDt = DvDt; - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - dvdt = obj.DvDt(obj.v); - [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); - state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); - end - - % Advances the solution vector one time step using the Runge-Kutta method given by - % obj.coeffs, updating the Residual Viscosity in each Runge-Kutta stage - function obj = step(obj) - obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.DvDt, obj.coeffs); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaRvBdf.m Fri Jun 28 13:13:17 2019 +0200 @@ -0,0 +1,93 @@ +classdef RungekuttaRvBdf < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + rkScheme % The particular RK scheme used for time integration + + + % Properties related to the residual viscositys + RV % Residual Viscosity operator + v_prev % Solution vector at previous time levels, used for the RV evaluation + DvDt % Function for computing the time deriative used for the RV evaluation + lowerBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation. + % dictates which accuracy the boot-strapping should start from. + upperBdfOrder % Orders of the approximation of the time deriative, used for the RV evaluation. + % Dictates the order of accuracy used once the boot-strapping is complete. + + + end + methods + function obj = RungekuttaRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + obj.RV = RV; + obj.lowerBdfOrder = bdfOrders.lowerBdfOrder; + obj.upperBdfOrder = bdfOrders.upperBdfOrder; + assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6)); + obj.v_prev = []; + obj.DvDt = rv.time.BdfDerivative(); + + if (rkOrder == 4) % Use specialized RK4 scheme + obj.rkScheme = @time.rk.rungekutta_4; + else + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(rkOrder); + coeffs = struct('s',s,'a',a,'b',b,'c',c); + obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); + end + + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + if (size(obj.v_prev,2) >= obj.lowerBdfOrder) + dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k); + [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); + else + viscosity = zeros(size(obj.v)); + dvdt = zeros(size(obj.v)); + Df = zeros(size(obj.v)); + firstOrderViscosity = zeros(size(obj.v)); + residualViscosity = zeros(size(obj.v)); + end + state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); + end + + function obj = step(obj) + nStoredStages = size(obj.v_prev,2); + + %Calculate viscosity for the new time level + if (nStoredStages >= obj.lowerBdfOrder) + viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k)); + else + viscosity = zeros(size(obj.v)); + end + + % Store current time level and update v_prev + if (nStoredStages < obj.upperBdfOrder) + obj.v_prev = [obj.v, obj.v_prev]; + else + obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); + obj.v_prev(:,1) = obj.v; + end + + % Fix the viscosity of the RHS function F + m = length(viscosity); + F_visc = @(v,t) obj.F(v, t, spdiags(viscosity,0,m,m)); + obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_visc); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaRvMultiGrid.m Fri Jun 28 13:13:17 2019 +0200 @@ -0,0 +1,65 @@ +classdef RungekuttaRvMultiGrid < time.Timestepper + properties + F % RHS of the ODE + F_coarse % RHS of the unstabalized ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + rkScheme % The particular RK scheme used for time integration + RV % Residual Viscosity operator + DvDt % Function for computing the time deriative used for the RV evaluation + v_unstable + viscosity + end + methods + + function obj = RungekuttaRvMultiGrid(F, F_coarse, k, t0, v0, RV, DvDt, order) + obj.F = F; + obj.F_coarse = F_coarse; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + + if (order == 4) % Use specialized RK4 scheme + obj.rkScheme = @time.rk.rungekutta_4; + else + % Extract the coefficients for the specified order + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(order); + coeffs = struct('s',s,'a',a,'b',b,'c',c); + obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs); + end + + obj.RV = RV; + obj.DvDt = DvDt; + obj.v_unstable = 0*v0; + obj.viscosity = 0*v0; + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + dvdt = obj.DvDt(obj.v_unstable); + [viscosity, Df, firstOrderViscosity, residualViscosity] = obj.RV.evaluate(obj.v, dvdt); + state = struct('v', obj.v, 'dvdt', dvdt, 'Df', Df, 'viscosity', obj.viscosity, 'residualViscosity', residualViscosity, 'firstOrderViscosity', firstOrderViscosity, 't', obj.t); + end + + % Advances the solution vector one time step using the Runge-Kutta method given by + % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps + function obj = step(obj) + m = length(obj.viscosity); + obj.v_unstable = obj.rkScheme(obj.v, obj.t, obj.k, obj.F_coarse); + obj.viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v_unstable)); + % Fix the viscosity of the stabilized RHS + F_stable = @(v,t) obj.F(v,t,spdiags(obj.viscosity,0,m,m)); + obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_stable); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file