Mercurial > repos > public > sbplib
changeset 1197:433c89bf19e0 feature/rv
Merge with default
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 07 Aug 2019 15:23:42 +0200 |
| parents | f6c571d8f22f (diff) 33c378e508d2 (current diff) |
| children | 5271c4670733 |
| files | +rv/+diffops/constructSymmetricD2.m +scheme/Beam2d.m +scheme/Burgers1d.m +scheme/Burgers2d.m +scheme/TODO.txt +scheme/Utux.m +scheme/Utux2d.m +scheme/Wave2dCurve.m +scheme/error1d.m +scheme/error2d.m +scheme/errorMax.m +scheme/errorRelative.m +scheme/errorSbp.m +scheme/errorVector.m +time/+cdiff/cdiff.m |
| diffstat | 32 files changed, 1540 insertions(+), 192 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/addClosuresToDiffOp.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,11 @@ +function D = addClosuresToDiffOp(D, closures) + for i = 1:size(closures,1) + for j = 1:size(closures,2) + closure = closures{i,j}; + if ~isa(closure, 'function_handle') + closure = @(v) closure*v; + end + D = @(v) D(v) + closure(v); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructDiffOps.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,136 @@ +function diffOpsStruct = constructDiffOps(method, varargin) + switch method + case {'standard'} + diffOpsStruct = diffOpsStandard(varargin{:}); + case {'bdf','backwards-difference-formula'} + diffOpsStruct = diffOpsBdf(varargin{:}); + case {'ms','multi-stage'} + diffOpsStruct = diffOpsMultiStage(varargin{:}); + case {'mg1','multi-grid1'} + diffOpsStruct = diffOpsMultiGrid1(varargin{:}); + case {'mg','mg2','multi-grid2'} % Default multigrid diffops + diffOpsStruct = diffOpsMultiGrid2(varargin{:}); + end +end + +function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) + % DiffOps for stabilized solution vector + [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + %% DiffOps for residual viscosity + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = @(v) -D_flux(v); + diffOpsStruct = struct('D_scheme', D_scheme,... + 'D_flux', D_flux,... + 'D_t', D_t,... + 'penalties_scheme', {penalties_scheme}); +end + +function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) + %% DiffOps for solution vector + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + %% DiffOps for residual viscosity + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = @(v) -D_flux(v); + diffOpsStruct = struct('D_scheme', D_scheme,... + 'D_flux', D_flux,... + 'penalties_scheme', {penalties_scheme}); +end + +% TODO: Remove? +function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) + % DiffOps for stabilized solution vector + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + % DiffOp for unstabilized solution vector + [D_unstable, closures] = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + D_unstable = rv.diffops.addClosuresToDiffOp(D_unstable, closures); + %% DiffOps for residual viscosity + D_t = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = @(v) -D_t(v); + % TODO: Use residual penalties as well here? + diffOpsStruct = struct('D_scheme', D_scheme,... + 'D_unstable', D_unstable,... + 'D_flux', D_flux,... + 'D_t', D_t,... + 'penalties_scheme', {penalties_scheme}); +end + +function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) + % DiffOps for stabilized solution vector + [D_scheme, penalties_scheme, D_f] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + % DiffOp for unstabilized solution vector + D_coarse = coarserGridDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); + %% DiffOps for residual viscosity + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs); + D_flux = @(v) -D_flux(v); + D_t = @(v) D_f(v); + % TODO: Use residual penalties as well here? + diffOpsStruct = struct('D_scheme', D_scheme,... + 'D_coarse', D_coarse,... + 'D_flux', D_flux,... + 'D_t', D_t,... + 'penalties_scheme', {penalties_scheme}); +end + +function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) + % DiffOps for stabilized solution vector + [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs); + % TODO: What orders to use here? + D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); + %% DiffOps for residual viscosity + D_flux = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder , schemeParams, opSet, BCs); + D_flux = @(v) -D_flux(v); + D_t = @(v) D_coarse(v); + diffOpsStruct = struct('D_scheme', D_scheme,... + 'D_flux', D_flux,... + 'D_t', D_t,... + 'penalties_scheme', {penalties_scheme}); +end + +% Note: Only works for equidistant grids. +function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) + m = g.m(); + lim = g.lim(); + m_c = (m-1)/2 + 1; + switch g.D() + case 1 + interpOps = sbp.InterpOpsOP(m_c, m, schemeOrder, schemeOrder); + I = interpOps.Iu2v.good; + g_c = grid.equidistant(m_c, lim{1}); + D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = @(v) I*D_c(v(1:2:end)); + case 2 + interpOps_x = sbp.InterpOpsOP(m_c(1), m(1), schemeOrder, schemeOrder); + I_x = interpOps_x.Iu2v.good; + interpOps_y = sbp.InterpOpsOP(m_c(2), m(2), schemeOrder, schemeOrder); + I_y = interpOps_y.Iu2v.good; + I = kron(I_x,I_y); + g_c = grid.equidistant(m_c, lim{1}, lim{2}); + ind = grid.funcToMatrix(g, 1:g.N()); + ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1); + D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = @(v) I*D_c(v(ind_c)); + end +end + +function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs) + m = g.m(); + lim = g.lim(); + m_c = (m-1)/2 + 1; + switch g.D() + case 1 + g_c = grid.equidistant(m_c, lim{1}); + [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = rv.diffops.addClosuresToDiffOp(D_c, closures); + x = g.x{1}; + x_c = x(1:2:end); + D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline'); + case 2 + g_c = grid.equidistant(m_c, lim{1}, lim{2}); + ind = grid.funcToMatrix(g, 1:g.N()); + ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1); + [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); + D_c = rv.diffops.addClosuresToDiffOp(D_c, closures); + D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(v(ind_c))),'spline')',g.N(),1); + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructFluxDiffOps.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,15 @@ +function [D_f, closures, penalties] = constructFluxDiffOps(scheme, g, order, schemeParams, opSet, BCs) + diffOp = scheme(g, order, schemeParams{:}, opSet); + if ~isa(diffOp.D, 'function_handle') + D_f = @(v) diffOp.D*v; + else + D_f = diffOp.D; + end + penalties = cell(size(BCs)); + closures = cell(size(BCs)); + for i = 1:size(BCs,1) + for j = 1:size(BCs,2) + [closures{i,j}, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructSchemeDiffOps.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,7 @@ +function [D_scheme, penalties, D_f] = constructSchemeDiffOps(scheme, g, order, schemeParams, opSet, BCs) + %% DiffOps for solution vector + [D_f, closures, penalties] = rv.diffops.constructFluxDiffOps(scheme, g, order, schemeParams, opSet, BCs); + D_scheme = rv.diffops.addClosuresToDiffOp(D_f, closures); + D2 = rv.diffops.constructSymmetricD2(g, order, opSet); + D_scheme = @(v,viscosity)(D_scheme(v) + D2(viscosity)*v); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+diffops/constructSymmetricD2.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,55 @@ +function D2 = constructSymmetricD2(g, order, opSet) + m = g.size(); + ops = cell(g.D(),1); + I = cell(g.D(),1); + for i = 1:g.D() + lim = {g.x{i}(1), g.x{i}(end)}; + ops{i} = opSet(m(i), lim, order); + I{i} = speye(m(i)); + end + + switch g.D() + case 1 + e_r = ops{1}.e_r; + e_l = ops{1}.e_l; + Hi = ops{1}.HI; + B = e_r*e_r' - e_l*e_l'; + if isequal(opSet,@sbp.D1Upwind) + Dm = ops{1}.Dm; + Dp = ops{1}.Dp; + M = Dm - Hi*B; + D2 = @(Viscosity) M*Viscosity*Dp; + else + % TODO: Fix Viscosity not being vector + d1_r = ops{1}.d1_r'; + d1_l = ops{1}.d1_l'; + D2 = @(Viscosity)ops{1}.D2(diag(Viscosity)) + Hi*(Viscosity(1,1)*e_l*d1_l - e_r*Viscosity(end,end)*d1_r); + end + case 2 + % TODO: + % Currently only implemented for upwind operators. + % Remove this part once the time-dependent D2 operator is implemented for other opSets + % or if it is decided that it should only be supported for upwind operators. + assert(isequal(opSet,@sbp.D1Upwind)) + e_e = kron(ops{1}.e_r,I{2}); + e_w = kron(ops{1}.e_l,I{2}); + Dm_x = kron(ops{1}.Dm,I{2}); + Dp_x = kron(ops{1}.Dp,I{2}); + H_x = kron(ops{1}.HI,I{2}); + B_x = e_e*e_e' - e_w*e_w'; + M_x = Dm_x-H_x*B_x; + D2_x = @(Viscosity) M_x*Viscosity*Dp_x; + + e_n = kron(I{1},ops{2}.e_r); + e_s = kron(I{1},ops{2}.e_l); + Dm_y = kron(I{1},ops{2}.Dm); + Dp_y = kron(I{1},ops{2}.Dp); + H_y = kron(I{1},ops{2}.HI); + B_y = e_n*e_n' - e_s*e_s'; + M_y = Dm_y-H_y*B_y; + D2_y = @(Viscosity) M_y*Viscosity*Dp_y; + D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity); + otherwise + error('3D not yet implemented') + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/BdfDerivative.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,39 @@ +classdef BdfDerivative < handle + properties + coefficients + end + methods + function obj = BdfDerivative() + obj.coefficients = obj.getBdfCoefficients(); + end + function DvDt = evaluate(obj, v, v_prev, dt) + order = size(v_prev,2); + DvDt = (obj.coefficients(order,1)*v - sum(obj.coefficients(order,2:order+1).*v_prev,2))/dt; + end + end + methods(Static) + function c = getBdfCoefficients() + c = zeros(9,10); + c(1,1) = 1; c(1,2) = 1; + c(2,1) = 3; c(2,2) = 4; c(2,3) = -1; + c(3,1) = 11; c(3,2) = 18; c(3,3) = -9; c(3,4) = 2; + c(4,1) = 25; c(4,2) = 48; c(4,3) = -36; c(4,4) = 16; c(4,5) = -3; + c(5,1) = 137; c(5,2) = 300; c(5,3) = -300; c(5,4) = 200; c(5,5) = -75; c(5,6) = 12; + c(6,1) = 147; c(6,2) = 360; c(6,3) = -450; c(6,4) = 400; c(6,5) = -225; c(6,6) = 72; c(6,7) = -10; + % Note: Higher orders than 6 are not stable, but can still be used as one-sided approximations of a derivatve + c(7,1) = 1089; c(7,2) = 2940; c(7,3) = -4410; c(7,4) = 4900; c(7,5) = -3675; c(7,6) = 1764; c(7,7) = -490; c(7,8) = 60; + c(8,1) = 2283; c(8,2) = 6720; c(8,3) = -11760; c(8,4) = 15680; c(8,5) = -14700; c(8,6) = 9408; c(8,7) = -3920; c(8,8) = 960; c(8,9) = -105; + c(9,1) = 7129; c(9,2) = 22680; c(9,3) = -45360; c(9,4) = 70560; c(9,5) = -79380; c(9,6) = 63504; c(9,7) = -35280; c(9,8) = 12960; c(9,9) = -2835; c(9,10) = 280; + + % Scale coefficients + c(2,:) = c(2,:)/2; + c(3,:) = c(3,:)/6; + c(4,:) = c(4,:)/12; + c(5,:) = c(5,:)/60; + c(6,:) = c(6,:)/60; + c(7,:) = c(7,:)/420; + c(8,:) = c(8,:)/840; + c(9,:) = c(9,:)/2520; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaRv.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,57 @@ +classdef RungekuttaRv < 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 = RungekuttaRv(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 stabilized RHS + m = length(obj.v); + F = @(v,t) obj.F(v,t,spdiags(obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v)),0,m,m)); + obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F); + 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 Wed Aug 07 15:23:42 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 <= 9)); + 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/RungekuttaRvInstage.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,79 @@ +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 + + % Takes one time step of size dt using the rungekutta method + % starting from v and where the function F(v,t,RV) gives the + % time derivatives. coeffs is a struct holding the RK coefficients + % for the specific method. RV is the residual viscosity which is updated + % in between the stages and after the updated solution is computed. + methods (Static) + function v = rungekuttaRV(v, t , dt, F, RV, DvDt, coeffs) + % Move one stage outside to avoid branching for updating the + % residual inside the loop. + k = zeros(length(v), coeffs.s); + k(:,1) = F(v,t,RV.evaluateViscosity(v,DvDt(v))); + + % Compute the intermediate stages k + for i = 2:coeffs.s + u = v; + for j = 1:i-1 + u = u + dt*coeffs.a(i,j)*k(:,j); + end + k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluateViscosity(u,DvDt(u))); + end + + % Compute the updated solution as a linear combination + % of the intermediate stages. + u = v; + for i = 1:coeffs.s + u = u + dt*coeffs.b(i)*k(:,i); + end + v = u; + end + + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaRvMultiGrid.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,68 @@ +classdef RungekuttaRvMultiGrid < time.Timestepper + properties + F % RHS of the ODE + F_coarse % RHS of the untabilized coarse grid 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_coarse + 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_coarse = 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_coarse); + [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) + % Fix the viscosity of the stabilized RHS + m = length(obj.viscosity); + F_stable = @(v,t) obj.F(v,t,spdiags(obj.viscosity,0,m,m)); + % Advance solution on unstabilized coarse mesh based on current solution + obj.v_coarse = obj.rkScheme(obj.v, obj.t, obj.k, obj.F_coarse); + % Advance solution on on stabilized mesh based on current viscosity + obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_stable); + % Compute viscosity for the next time time level using the advanced solution + obj.viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v_coarse)); + 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/RungekuttaRvMultiStage.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,70 @@ +classdef RungekuttaRvMultiStage < time.Timestepper + properties + F % RHS of the ODE + F_unstable % RHS of the unstabilized 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) + + % Advance solution by unstabilized scheme + obj.v_unstable = obj.rkScheme(obj.v, obj.t, obj.k, obj.F_unstable); + % Compute viscosity for current time level based on unstable solution (from next time level) + % and the current solution. + obj.viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v_unstable)); + % Fix the viscosity of the stabilized RHS + m = length(obj.viscosity); + F_stable = @(v,t) obj.F(v,t,spdiags(obj.viscosity,0,m,m)); + % Advance solutiont to next time level by stabilized scheme. + 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/getRvTimestepper.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,70 @@ +function ts = getRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + switch opt.method + case 'rkRv' + ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0); + case 'rkRvBdf' + ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0); + case 'rkRvMs' + ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0); + case 'rkRvMg' + ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0); + case 'rkRvInstage' + ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0); + otherwise + error('Timestepping method ''%s'' not supported',method); + end +end + +function fhData = dataFunctionHandle(data) + if isa(data, 'function_handle') + switch nargin(data) + case 1 + fhData = @(v, t) data(v); + case 2 + fhData = @(v, t) data(v,t); + otherwise + error('Incorrect number of input arguments'); + end + else + fhData = @(v, t) data; + end +end + +function F = stabilizedRhs(D, data) + fhData = dataFunctionHandle(data); + F = @(v, t, viscosity) D(v, viscosity) + fhData(v,t); +end + +function F = unstabilizedRhs(D, data) + fhData = dataFunctionHandle(data); + F = @(v, t) D(v) + fhData(v,t); +end + +function ts = rkRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + F = stabilizedRhs(diffOpStruct.D_scheme, data); + ts = rv.time.RungekuttaRv(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder); +end + +function ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + F = stabilizedRhs(diffOpStruct.D_scheme, data); + ts = rv.time.RungekuttaRvBdf(F, opt.k, t0, v0, residualViscosity, opt.rkOrder, opt.bdfOrders); +end + +function ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + F = stabilizedRhs(diffOpStruct.D_scheme, data); + F_unstab = unstabilizedRhs(diffOpStruct.D_unstable, data); + ts = rv.time.RungekuttaRvMultiStage(F, F_unstab, opt.k, t0, v0,... + residualViscosity, diffOpStruct.D_t, opt.rkOrder); +end + +function ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + F = stabilizedRhs(diffOpStruct.D_scheme, data); + F_coarse = unstabilizedRhs(diffOpStruct.D_coarse, data); + ts = rv.time.RungekuttaRvMultiGrid(F, F_coarse, opt.k, t0, v0,... + residualViscosity, diffOpStruct.D_t, opt.rkOrder); +end + +function ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0) + F = stabilizedRhs(diffOpStruct.D_scheme, data); + ts = rv.time.RungekuttaRvInstage(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/ResidualViscosity.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,110 @@ +classdef ResidualViscosity < handle + properties + g % grid + Df % Diff op approximating the gradient of the flux f(u) + waveSpeed % Wave speed at each grid point, e.g f'(u). %TBD: Better naming? + Cmax % Constant controlling relative amount of upwind dissipation + Cres % Constant controling relative amount of upwind dissipation + h % Length scale used for scaling the viscosity. Typically grid spacing. + normalization % Function used to normalize the residual such that it is amplified in the + % shocks and suppressed elsewhere. + Mres % Coefficients for the residual viscosity + Mfirst % Coefficients for the first order viscosity + fRes % Function handle for computing the residual. + end + + methods + % TODO: pass opt struct with waveSpeed, normalization etc. + % TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without + % wrapping it in a function. + function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization, postProcess) + %default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf))); + obj.Df = Df; + obj.waveSpeed = waveSpeed; + obj.h = h; + obj.Cmax = Cmax; + + obj.Cres = Cres; + obj.normalization = normalization; + obj.g = g; + obj.Mres = obj.Cres*obj.h^2; + obj.Mfirst = obj.Cmax*obj.h; + switch postProcess + case {'', 'none'} + obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); + % TBD: Keep? + % case {'filt', 'filter'} + % order = 4; + % F = obj.shapiroFilter(obj.g, order); + % obj.Mres = F*obj.Mres; + % obj.fRes = @(v,dvdt) obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); + case {'max', 'maximum neighbors'} + switch g.D() + case 1 + obj.fRes = @(v,dvdt) movmax(obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v),3); + case 2 + obj.fRes = @obj.maxResidualNeighbors2d; + end + end + end + + function viscosity = evaluateViscosity(obj, v, dvdt) + viscosity = min(obj.Mfirst*abs(obj.waveSpeed(v)), obj.fRes(v,dvdt)); + end + + function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) + Df = obj.Df(v); + firstOrderViscosity = obj.Mfirst*abs(obj.waveSpeed(v)); + residualViscosity = obj.fRes(v,dvdt); + viscosity = min(firstOrderViscosity, residualViscosity); + end + + function res = maxResidualNeighbors2d(obj,v,dvdt) + res = obj.Mres*abs(dvdt + obj.Df(v))./obj.normalization(v); + resMat = grid.funcToMatrix(obj.g,res); + res = reshape(max(movmax(resMat,3,1),movmax(resMat,3,2))',obj.g.N(),1); + end + end + methods (Static) + function minmaxDiff = minmaxDiffNeighborhood1d(u) + umax = movmax(u,3); + umin = movmin(u,3); + minmaxDiff = umax - umin; + end + + function minmaxDiff = minmaxDiffNeighborhood2d(g, u) + uMatrix = grid.funcToMatrix(g,u); + umax = max(movmax(uMatrix,3,1),movmax(uMatrix,3,2)); + umin = min(movmin(uMatrix,3,1),movmin(uMatrix,3,2)); + minmaxDiff = umax - umin; + minmaxDiff = reshape(minmaxDiff',g.N(),1); + end + function F = shapiroFilter(g, order) + switch order + case 2 + F = spdiags(repmat([1 2 1],g.m(1),1), -1:1, g.m(1), g.m(1)); + case 4 + F = spdiags(repmat([-1 4 10 4 -1],g.m(1),1), -2:2, g.m(1), g.m(1)); + case 6 + F = spdiags(repmat([1 -6 15 44 15 -6 1],g.m(1),1), -3:3, g.m(1), g.m(1)); + case 8 + F = spdiags(repmat([-1 8 -28 56 186 56 -28 8 -1],g.m(1),1), -4:4, g.m(1), g.m(1)); + end + F = 1/2^order * F; + % Fx = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(1),1), -1:1, g.m(1), g.m(1)); + % Fx(1,:) = 0; + % Fx(1,1:2) = 1/(mid+1)*[mid 1]; + % Fx(end,:) = 0; + % Fx(end,end-1:end) = 1/(mid+1)*[1 mid]; + % Fy = spdiags(repmat(1/(mid+2)*[1 mid 1],g.m(2),1), -1:1, g.m(2), g.m(2)); + % Fy(1,:) = 0; + % Fy(1,1:2) = 1/(mid+1)*[mid 1]; + % Fy(end,:) = 0; + % Fy(end,end-1:end) = 1/(mid+1)*[1 mid]; + % Ix = speye(g.m(1)); + % Iy = speye(g.m(1)); + % F = kron(Fx, Iy) + kron(Ix, Fy); + end + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/dissipationOperator.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,67 @@ +%% Function that constructs artificial dissipation operators using undivided differences +function D = dissipationOperator(m, order, Hinv, scaling) + % TBD: Add or remove D_2 and Dp/Dm? + % d2=[1 2 1]; + % D_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ... + % diag(ones(m-1,1),1)); + % D_2(1,1:3)=[d2];D_2(m,m-2:m)=[d2]; + % %Dm + % DD_m=(diag(ones(m-1,1),+1)-diag(ones(m,1),0)); + % DD_m(m,m-1:m)=[-1 1]; + % %Dp + % DD_p=(-diag(ones(m-1,1),-1)+diag(ones(m,1),0)); + % DD_p(1,1:2)=[-1 1]; + + switch order + case 1 + DD_1=(diag(ones(m-1,1),+1)-diag(ones(m,1),0)); + DD_1(m,m-1:m)=[-1 1]; + D = DD_2'*DD_2; + case 2 + dd2=0*[1 -2 1]; + DD_2=(diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)+ ... + diag(ones(m-1,1),1)); + DD_2(1,1:3)=[dd2];DD_2(m,m-2:m)=[dd2]; + D = DD_2'*DD_2; + case 3 + d3=0*[-1 3 -3 1]; + DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ... + diag(ones(m-1,1),1)); + DD_3(1:2,1:4)=[d3;d3]; + DD_3(m,m-3:m)=[d3]; + D = DD_3'*DD_3; + case 4 + default_arg('scaling', 1/12); + d4=0*[1 -4 6 -4 1]; + DD_4=(diag(ones(m-2,1),2)-4*diag(ones(m-1,1),1)+6*diag(ones(m,1),0)-4*diag(ones(m-1,1),-1)+diag(ones(m-2,1),-2)); + DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4]; + D = DD_4'*DD_4; + case 5 + d5=0*[-1 5 -10 10 -5 1]; + DD_5=(-diag(ones(m-3,1),-3)+5*diag(ones(m-2,1),-2)-10*diag(ones(m-1,1),-1)+10*diag(ones(m,1),0)-5*diag(ones(m-1,1),1)+diag(ones(m-2,1),2)); + DD_5(1:3,1:6)=[d5;d5;d5]; + DD_5(m-1:m,m-5:m)=[d5;d5]; + D = DD_5'*DD_5; + case 6 + default_arg('scaling', 1/60); + d6=0*[1 -6 15 -20 15 -6 1]; + DD_6=(diag(ones(m-3,1),3)-6*diag(ones(m-2,1),2)+15*diag(ones(m-1,1),1)-20*diag(ones(m,1),0)+15*diag(ones(m-1,1),-1)-6*diag(ones(m-2,1),-2)+diag(ones(m-3,1),-3)); + DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6]; + D = DD_6'*DD_6; + case 7 + d7=0*[-1 7 -21 35 -35 21 -7 1]; + DD_7=(-diag(ones(m-4,1),-4)+7*diag(ones(m-3,1),-3)-21*diag(ones(m-2,1),-2)+35*diag(ones(m-1,1),-1)-35*diag(ones(m,1),0)+21*diag(ones(m-1,1),1)-7*diag(ones(m-2,1),2)+diag(ones(m-3,1),3)); + DD_7(1:4,1:8)=[d7;d7;d7;d7]; + DD_7(m-2:m,m-7:m)=[d7;d7;d7]; + D = DD_7'*DD_7; + case 9 + d9=0*[-1 9 -36 84 -126 126 -84 36 -9 1]; + DD_9=(-diag(ones(m-5,1),-5)+9*diag(ones(m-4,1),-4)-36*diag(ones(m-3,1),-3)+84*diag(ones(m-2,1),-2)-126*diag(ones(m-1,1),-1)+126*diag(ones(m,1),0)-84*diag(ones(m-1,1),1)+36*diag(ones(m-2,1),2)-9*diag(ones(m-3,1),3)+diag(ones(m-4,1),4)); + DD_9(1:5,1:10)=[d9;d9;d9;d9;d9]; + DD_9(m-3:m,m-9:m)=[d9;d9;d9;d9]; + D = DD_9'*DD_9; + otherwise + error('Order not yet supported', order); + end + D = scaling*sparse(Hinv*D); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Burgers1d.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,145 @@ +classdef Burgers1d < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + grid % Grid + order % Order accuracy for the approximation + + H % Discrete norm + D + + D1 + Hi + e_l + e_r + end + + methods + function obj = Burgers1d(g, order, pde_form, fluxSplitting, opSet) + default_arg('opSet',@sbp.D2Standard); + default_arg('fluxSplitting',@(v)max(abs(v))); + assertType(g, 'grid.Cartesian'); + + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + + ops = opSet(m, xlim, order); + + if (isequal(opSet, @sbp.D1Upwind)) + obj.D1 = (ops.Dp + ops.Dm)/2; + DissOp = (ops.Dm - ops.Dp)/2; + switch pde_form + case 'quasi-linear' + obj.D = @(v) -((spdiag(v)*obj.D1 + fluxSplitting(v)*DissOp)*v); + case 'skew-symmetric' + obj.D = @(v) -(1/3*obj.D1*(v.*v) + (1/3*spdiag(v)*obj.D1 + fluxSplitting(v)*DissOp)*v); + case 'conservative' + obj.D = @(v) -(1/2*obj.D1*(v.*v) + fluxSplitting(v)*DissOp*v); + end + else + obj.D1 = ops.D1; + switch pde_form + case 'quasi-linear' + obj.D = @(v) -(spdiag(v)*obj.D1*v); + case 'skew-symmetric' + obj.D = @(v) -(1/3*obj.D1*(v.*v) + 1/3*spdiag(v)*obj.D1*v); + case 'conservative' + obj.D = @(v) -1/2*obj.D1*(v.*v); + end + end + obj.grid = g; + + obj.H = ops.H; + obj.Hi = ops.HI; + + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; + + obj.m = m; + obj.h = ops.h; + obj.order = order; + end + + % Closure functions return the operators applied to the own doamin to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + function [closure, penalty] = boundary_condition(obj, boundary, type) + default_arg('type','dirichlet'); + s = obj.getBoundarySign(boundary); + e = obj.getBoundaryOperator('e', boundary); + index = obj.getBoundaryIndex(boundary); + switch type + % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3 + % with +- at left/right boundaries + case {'D', 'd', 'dirichlet', 'Dirichlet'} + % tau = s*e; + % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index)); + % penalty = -obj.Hi*tau; + + penalty_parameter = 1/3; + tau = @(v) s*penalty_parameter*obj.Hi*e*(v(index)-s*abs(v(index)))/2; + closure = @(v) tau(v)*v(index); + penalty = @(v) -tau(v); + otherwise + error('No such boundary condition: type = %s',type); + end + end + + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + + switch boundary + case {'r'} + s = 1; + case {'l'} + s = -1; + end + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + % boundary -- string + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e'}) + assertIsMember(boundary, {'l', 'r'}) + + o = obj.([op, '_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + H_b = 1; + end + + % Returns the boundary index. The right boundary has the last index + % boundary -- string + function index = getBoundaryIndex(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + switch boundary + case {'r'} + index = length(obj.e_r); + case {'l'} + index = 1; + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.grid.m; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Burgers2d.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,194 @@ +classdef Burgers2d < scheme.Scheme + properties + grid % Physical grid + order % Order accuracy for the approximation + + D % Non-stabilized scheme operator + H % Discrete norm + H_x, H_y % Norms in the x and y directions + Hi % Kroneckered norm inverse + % Boundary operators + e_w, e_e, e_s, e_n + end + + methods + function obj = Burgers2d(g, order, pde_form, fluxSplitting, opSet) + default_arg('opSet',@sbp.D2Standard); + default_arg('fluxSplitting',{@(v)max(abs(v)),@(v)max(abs(v))}); + assertType(g, 'grid.Cartesian'); + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + obj.grid = g; + + % Operator sets + ops_x = opSet(m_x, xlim, order); + ops_y = opSet(m_y, ylim, order); + Ix = speye(m_x); + Iy = speye(m_y); + + % Norms + Hx = ops_x.H; + Hy = ops_y.H; + Hxi = ops_x.HI; + Hyi = ops_y.HI; + + obj.H_x = Hx; + obj.H_y = Hy; + obj.H = kron(Hx,Hy); + obj.Hi = kron(Hxi,Hyi); + + % Derivatives + if (isequal(opSet,@sbp.D1Upwind)) + Dx = kron((ops_x.Dp + ops_x.Dm)/2,Iy); + Dy = kron(Ix,(ops_y.Dp + ops_y.Dm)/2); + DissOpx = kron((ops_x.Dp - ops_x.Dm)/2,Iy); + DissOpy = kron(Ix,(ops_y.Dp - ops_y.Dm)/2); + D1 = Dx + Dy; + switch pde_form + case 'skew-symmetric' + D = -1/3*D1; + switch length(fluxSplitting) + case 1 + DissOp = DissOpx + DissOpy; + obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOp)*v; + case 2 + obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v; + end + case 'conservative' + D = -1/2*D1; + switch length(fluxSplitting) + case 1 + DissOp = DissOpx + DissOpy; + % TODO: Check if we can use fluxSplitting{1} here instead + obj.D = @(v) D*(v.*v) + max(abs(v))*DissOp*v; + case 2 + obj.D = @(v) D*(v.*v) + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v; + end + + end + else + Dx = kron(ops_x.D1,Iy); + Dy = kron(Ix,ops_y.D1); + D1 = Dx + Dy; + switch pde_form + case 'skew-symmetric' + D = -1/3*D1; + obj.D = @(v) D*(v.*v) + spdiags(v,0,m_tot,m_tot)*D*v; + case 'conservative' + D = -1/2*D1; + obj.D = @(v) D*(v.*v); + end + end + + % Boundary operators + obj.e_w = kr(ops_x.e_l, Iy); + obj.e_e = kr(ops_x.e_r, Iy); + obj.e_s = kr(Ix, ops_y.e_l); + obj.e_n = kr(Ix, ops_y.e_r); + + obj.order = order; + end + + % Closure functions return the operators applied to the own doamin to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); + s = obj.getBoundarySign(boundary); + e = obj.getBoundaryOperator('e', boundary); + indices = obj.getBoundaryIndices(boundary); + H_1d = obj.getOneDirectionalNorm(boundary); + switch type + case {'D', 'd', 'dirichlet', 'Dirichlet'} + penalty_parameter = 1/3; + Tau = s*penalty_parameter*obj.Hi*e*H_1d/2; + m = obj.grid.m; + tau = @(v) Tau*spdiags((v(indices)-s*abs(v(indices))),0,m(1),m(2)); + closure = @(v) Tau*((v(indices)-s*abs(v(indices))).*v(indices)); + penalty = @(v) -tau(v); + otherwise + error('No such boundary condition: type = %s',type); + end + + + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} + s = -1; + end + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + % boundary -- string + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e'}) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + o = obj.([op, '_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + H_b = obj.(['H_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_1d = getOneDirectionalNorm(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary + case {'w','e'} + H_1d = obj.H_y; + case {'s','n'} + H_1d = obj.H_x; + end + end + + % Returns the indices of the boundary points in the grid matrix + % boundary -- string + function I = getBoundaryIndices(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.grid.m)); + switch boundary + case 'w' + I = ind(1,:); + case 'e' + I = ind(end,:); + case 's' + I = ind(:,1)'; + case 'n' + I = ind(:,end)'; + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.grid.m; + end + end +end \ No newline at end of file
--- a/+scheme/Utux.m Thu May 02 13:25:14 2019 -0700 +++ b/+scheme/Utux.m Wed Aug 07 15:23:42 2019 +0200 @@ -5,6 +5,9 @@ grid % Grid order % Order accuracy for the approximation + a % Wave speed + % Can either be a constant or function handle. + H % Discrete norm D @@ -12,13 +15,21 @@ Hi e_l e_r - v0 end methods - function obj = Utux(g, order, opSet) + function obj = Utux(g, order, a, fluxSplitting, opSet) default_arg('opSet',@sbp.D2Standard); + default_arg('a',1); + default_arg('fluxSplitting',[]); + + assertType(g, 'grid.Cartesian'); + if isa(a, 'function_handle') + obj.a = spdiag(grid.evalOn(g, a)); + else + obj.a = a; + end m = g.size(); xl = g.getBoundary('l'); @@ -26,7 +37,15 @@ xlim = {xl, xr}; ops = opSet(m, xlim, order); - obj.D1 = ops.D1; + + if (isequal(opSet, @sbp.D1Upwind)) + obj.D1 = (ops.Dp + ops.Dm)/2; + DissOp = (ops.Dm - ops.Dp)/2; + obj.D = -(obj.a*obj.D1 + fluxSplitting*DissOp); + else + obj.D1 = ops.D1; + obj.D = -obj.a*obj.D1; + end obj.grid = g; @@ -35,12 +54,10 @@ obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - end % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. @@ -51,27 +68,53 @@ % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dirichlet'); - tau =-1*obj.e_l; - closure = obj.Hi*tau*obj.e_l'; + s = obj.getBoundarySign(boundary); + e = obj.getBoundaryOperator('e', boundary); + switch boundary + % Can only specify boundary condition where there is inflow + % Extract the postivie resp. negative part of a, for the left + % resp. right boundary, and set other values of a to zero. + % Then the closure will effectively only contribute to inflow boundaries + case {'l'} + a_inflow = obj.a; + a_inflow(a_inflow < 0) = 0; + case {'r'} + a_inflow = obj.a; + a_inflow(a_inflow > 0) = 0; + end + tau = s*a_inflow*e; + closure = obj.Hi*tau*e'; penalty = -obj.Hi*tau; - end function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) switch boundary % Upwind coupling case {'l','left'} - tau = -1*obj.e_l; + tau = -1*obj.a*obj.e_l; closure = obj.Hi*tau*obj.e_l'; penalty = -obj.Hi*tau*neighbour_scheme.e_r'; case {'r','right'} - tau = 0*obj.e_r; + tau = 0*obj.a*obj.e_r; closure = obj.Hi*tau*obj.e_r'; penalty = -obj.Hi*tau*neighbour_scheme.e_l'; end end + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + + switch boundary + case {'r'} + s = 1; + case {'l'} + s = -1; + end + end + % Returns the boundary operator op for the boundary specified by the string boundary. % op -- string % boundary -- string
--- a/+scheme/Utux2d.m Thu May 02 13:25:14 2019 -0700 +++ b/+scheme/Utux2d.m Wed Aug 07 15:23:42 2019 +0200 @@ -4,7 +4,6 @@ h % Grid spacing grid % Grid order % Order accuracy for the approximation - v0 % Initial data a % Wave speed a = [a1, a2]; % Can either be a constant vector or a cell array of function handles. @@ -25,10 +24,11 @@ methods - function obj = Utux2d(g ,order, opSet, a) + function obj = Utux2d(g ,order, a, fluxSplitting, opSet) default_arg('a',1/sqrt(2)*[1, 1]); default_arg('opSet',@sbp.D2Standard); + default_arg('fluxSplitting',[]); assertType(g, 'grid.Cartesian'); if iscell(a) @@ -74,10 +74,25 @@ obj.Hyi = kron(Ix,Hyi); % Derivatives - Dx = ops_x.D1; - Dy = ops_y.D1; - obj.Dx = kron(Dx,Iy); - obj.Dy = kron(Ix,Dy); + if (isequal(opSet,@sbp.D1Upwind)) + Dx = (ops_x.Dp + ops_x.Dm)/2; + Dy = (ops_y.Dp + ops_y.Dm)/2; + obj.Dx = kron(Dx,Iy); + obj.Dy = kron(Ix,Dy); + DissOpx = (ops_x.Dm - ops_x.Dp)/2; + DissOpy = (ops_y.Dm - ops_y.Dp)/2; + DissOpx = kron(DissOpx,Iy); + DissOpy = kron(Ix,DissOpy); + + obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*DissOpx + fluxSplitting{2}*DissOpy); + else + Dx = ops_x.D1; + Dy = ops_y.D1; + obj.Dx = kron(Dx,Iy); + obj.Dy = kron(Ix,Dy); + + obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); + end % Boundary operators obj.e_w = kr(ops_x.e_l, Iy); @@ -89,31 +104,40 @@ obj.h = [ops_x.h ops_y.h]; obj.order = order; obj.a = a; - obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); - end % Closure functions return the opertors applied to the own domain to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. + % type is a string specifying the type of boundary condition if there are several. %TBD Remove type here? Only dirichlet applicable? % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dirichlet'); - - sigma = -1; % Scalar penalty parameter + s = obj.getBoundarySign(boundary); + e = obj.getBoundaryOperator('e', boundary); + H_1d = obj.getOneDirectionalNorm(boundary); switch boundary + % Can only specify boundary condition where there is inflow + % Extract the postivie resp. negative part of a, for the left + % resp. right boundaries, and set other values of a to zero. + % Then the closure will effectively only contribute to inflow boundaries case {'w','W','west','West'} - tau = sigma*obj.a{1}*obj.e_w*obj.H_y; - closure = obj.Hi*tau*obj.e_w'; - + a_inflow = obj.a{1}; + a_inflow(a_inflow < 0) = 0; + case {'e','E','east','East'} + a_inflow = obj.a{1}; + a_inflow(a_inflow > 0) = 0; case {'s','S','south','South'} - tau = sigma*obj.a{2}*obj.e_s*obj.H_x; - closure = obj.Hi*tau*obj.e_s'; + a_inflow = obj.a{2}; + a_inflow(a_inflow < 0) = 0; + case {'n','N','north','North'} + a_inflow = obj.a{2}; + a_inflow(a_inflow > 0) = 0; end + tau = s*a_inflow*e*H_1d; + closure = obj.Hi*tau*e'; penalty = -obj.Hi*tau; - end % type Struct that specifies the interface coupling. @@ -277,6 +301,18 @@ end + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} + s = -1; + end + end + % Returns the boundary operator op for the boundary specified by the string boundary. % op -- string % boundary -- string @@ -297,6 +333,20 @@ H_b = obj.(['H_', boundary]); end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_1d = getOneDirectionalNorm(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary + case {'w','e'} + H_1d = obj.H_y; + case {'s','n'} + H_1d = obj.H_x; + end + end + function N = size(obj) N = obj.m; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/butcherTableau.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,46 @@ +function [s,a,b,c] = butcherTableau(order) +% TODO: Change order from a double to string. +switch order + + case 3 + % TVD (Total Variational Diminishing) + s = 3; + a = zeros(s,s-1); + a(2,1) = 1; + a(3,1) = 1/4; a(3,2) = 1/4; + b = [1/6, 1/6, 2/3]; + c = [0 1 1/2]; + case 4 + % Standard RK4 + s = 4; + a = zeros(s,s-1); + a(2,1) = 1/2; + a(3,1) = 0; a(3,2) = 1/2; + a(4,1) = 0; a(4,2) = 0; a(4,3) = 1; + b = [1/6 1/3 1/3 1/6]; + c = [0, 1/2, 1/2, 1]; + % case 4-3/8 + % % 3/8 RK4 (Kuttas method). Lower truncation error, more flops + % s = 4; + % a = zeros(s,s-1); + % a(2,1) = 1/3; + % a(3,1) = -1/3; a(3,2) = 1; + % a(4,1) = 1; a(4,2) = -1; a(4,3) = 1; + % b = [1/8 3/8 3/8 1/8]; + % c = [0, 1/3, 2/3, 1]; + case 6 + % Runge-Kutta 6 from Alshina07 + s = 7; + a = zeros(s,s-1); + a(2,1) = 4/7; + a(3,1) = 115/112; a(3,2) = -5/16; + a(4,1) = 589/630; a(4,2) = 5/18; a(4,3) = -16/45; + a(5,1) = 229/1200 - 29/6000*sqrt(5); a(5,2) = 119/240 - 187/1200*sqrt(5); a(5,3) = -14/75 + 34/375*sqrt(5); a(5,4) = -3/100*sqrt(5); + a(6,1) = 71/2400 - 587/12000*sqrt(5); a(6,2) = 187/480 - 391/2400*sqrt(5); a(6,3) = -38/75 + 26/375*sqrt(5); a(6,4) = 27/80 - 3/400*sqrt(5); a(6,5) = (1+sqrt(5))/4; + a(7,1) = -49/480 + 43/160*sqrt(5); a(7,2) = -425/96 + 51/32*sqrt(5); a(7,3) = 52/15 - 4/5*sqrt(5); a(7,4) = -27/16 + 3/16*sqrt(5); a(7,5) = 5/4 - 3/4*sqrt(5); a(7,6) = 5/2 - 1/2*sqrt(5); + b = [1/12 0 0 0 5/12 5/12 1/12]; + c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; + otherwise + error('That Runge-Kutta order is not implemented', order) + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/get_rk4_time_step.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,21 @@ +% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values. +function k = get_rk4_time_step(lambda,l_type) + default_arg('l_type','complex') + + rad = abs(lambda); + if strcmp(l_type,'real') + % Real eigenvalue + % kl > -2.7852 + k = 2.7852/rad; + + elseif strcmp(l_type,'imag') + % Imaginary eigenvalue + % |kl| < 2.8284 + k = 2.8284/rad; + elseif strcmp(l_type,'complex') + % |kl| < 2.5 + k = 2.5/rad; + else + error('l_type must be one of ''real'',''imag'' or ''complex''.') + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rk4_stability.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,58 @@ +function rk_stability() + ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); + circ = @(z)(abs(z)); + + + % contour(X,Y,z) + ax = [-4 2 -3 3]; + % hold on + fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2]) + hold on + r = 2.6; + fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r') + hold off + % contour(X,Y,z,[1,1],'b') + axis(ax) + title('4th order Runge-Kutta stability region') + xlabel('Re') + ylabel('Im') + axis equal + grid on + box on + hold off + % surf(X,Y,z) + + + rk4roots() +end + +function fcontour(f,levels,x_lim,y_lim,opt) + default_arg('opt','b') + x = linspace(x_lim(1),x_lim(2)); + y = linspace(y_lim(1),y_lim(2)); + [X,Y] = meshgrid(x,y); + mu = X+ 1i*Y; + + z = f(mu); + + contour(X,Y,z,levels,opt) + +end + + +function rk4roots() + ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); + % Roots for real evalues: + F = @(x)(abs(ruku4(x))-1); + real_x = fzero(F,-3); + + % Roots for imaginary evalues: + F = @(x)(abs(ruku4(1i*x))-1); + imag_x1 = fzero(F,-3); + imag_x2 = fzero(F,3); + + + fprintf('Real x = %f\n',real_x) + fprintf('Imag x = %f\n',imag_x1) + fprintf('Imag x = %f\n',imag_x2) +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,20 @@ +% Takes one time step of size dt using the rungekutta method +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. coeffs is a struct holding the RK coefficients +% for the specific method. +function v = rungekutta(v, t , dt, F, coeffs) + % Compute the intermediate stages k + k = zeros(length(v), coeffs.s); + for i = 1:coeffs.s + u = v; + for j = 1:i-1 + u = u + dt*coeffs.a(i,j)*k(:,j); + end + k(:,i) = F(u,t+coeffs.c(i)*dt); + end + % Compute the updated solution as a linear combination + % of the intermediate stages. + for i = 1:coeffs.s + v = v + dt*coeffs.b(i)*k(:,i); + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/rungekutta_4.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,10 @@ +% Takes one time step of size dt using the rungekutta method +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = rungekutta_4(v, t , dt, F) + k1 = F(v ,t ); + k2 = F(v+0.5*dt*k1,t+0.5*dt); + k3 = F(v+0.5*dt*k2,t+0.5*dt); + k4 = F(v+ dt*k3,t+ dt); + v = v + (1/6)*(k1+2*(k2+k3)+k4)*dt; +end \ No newline at end of file
--- a/+time/+rk4/get_rk4_time_step.m Thu May 02 13:25:14 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,21 +0,0 @@ -% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values. -function k = get_rk4_time_step(lambda,l_type) - default_arg('l_type','complex') - - rad = abs(lambda); - if strcmp(l_type,'real') - % Real eigenvalue - % kl > -2.7852 - k = 2.7852/rad; - - elseif strcmp(l_type,'imag') - % Imaginary eigenvalue - % |kl| < 2.8284 - k = 2.8284/rad; - elseif strcmp(l_type,'complex') - % |kl| < 2.5 - k = 2.5/rad; - else - error('l_type must be one of ''real'',''imag'' or ''complex''.') - end -end \ No newline at end of file
--- a/+time/+rk4/rk4_stability.m Thu May 02 13:25:14 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,58 +0,0 @@ -function rk_stability() - ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); - circ = @(z)(abs(z)); - - - % contour(X,Y,z) - ax = [-4 2 -3 3]; - % hold on - fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2]) - hold on - r = 2.6; - fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r') - hold off - % contour(X,Y,z,[1,1],'b') - axis(ax) - title('4th order Runge-Kutta stability region') - xlabel('Re') - ylabel('Im') - axis equal - grid on - box on - hold off - % surf(X,Y,z) - - - rk4roots() -end - -function fcontour(f,levels,x_lim,y_lim,opt) - default_arg('opt','b') - x = linspace(x_lim(1),x_lim(2)); - y = linspace(y_lim(1),y_lim(2)); - [X,Y] = meshgrid(x,y); - mu = X+ 1i*Y; - - z = f(mu); - - contour(X,Y,z,levels,opt) - -end - - -function rk4roots() - ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4)); - % Roots for real evalues: - F = @(x)(abs(ruku4(x))-1); - real_x = fzero(F,-3); - - % Roots for imaginary evalues: - F = @(x)(abs(ruku4(1i*x))-1); - imag_x1 = fzero(F,-3); - imag_x2 = fzero(F,3); - - - fprintf('Real x = %f\n',real_x) - fprintf('Imag x = %f\n',imag_x1) - fprintf('Imag x = %f\n',imag_x2) -end \ No newline at end of file
--- a/+time/+rk4/rungekutta_4.m Thu May 02 13:25:14 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -% 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_4(v, t , k, F) - k1 = F(v ,t ); - k2 = F(v+0.5*k*k1,t+0.5*k); - k3 = F(v+0.5*k*k2,t+0.5*k); - k4 = F(v+ k*k3,t+ k); - v = v + (1/6)*(k1+2*(k2+k3)+k4)*k; -end \ No newline at end of file
--- a/+time/+rk4/rungekutta_6.m Thu May 02 13:25:14 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -% 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_6(v, t , k, F) - s = 7 - k = zeros(length(v),s) - a = zeros(7,6); - c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; - b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12]; - a = [ - 0, 0, 0, 0, 0, 0; - 4/7, 0, 0, 0, 0, 0; - 115/112, -5/16, 0, 0, 0, 0; - 589/630, 5/18, -16/45, 0, 0, 0; - 229/1200 - 29/6000*sqrt(5), 119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5), 0, 0; - 71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4, 0; - -49/480 + 43/160*sqrt(5), -425/96 + 51/32*sqrt(5), 52/15 - 4/5*sqrt(5), -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5); - ] - - for i = 1:s - u = v - for j = 1: i-1 - u = u + h*a(i,j) * k(:,j) - end - k(:,i) = F(t+c(i)*k,u) - end - - for i = 1:s - v = v + k*b(i)*k(:,i) - end -end
--- a/+time/Rk4SecondOrderNonlin.m Thu May 02 13:25:14 2019 -0700 +++ b/+time/Rk4SecondOrderNonlin.m Wed Aug 07 15:23:42 2019 +0200 @@ -61,7 +61,7 @@ end function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Rungekutta.m Wed Aug 07 15:23:42 2019 +0200 @@ -0,0 +1,46 @@ +classdef Rungekutta < time.Timestepper + properties + F % RHS of the ODE + k % Time step + t % Time point + v % Solution vector + n % Time level + scheme % The scheme used for the time stepping, e.g rk4, rk6 etc. + end + + + methods + % Timesteps v_t = F(v,t), using RK with specfied order from t = t0 with + % timestep k and initial conditions v = v0 + function obj = Rungekutta(F, k, t0, v0, order) + default_arg('order',4); + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % TBD: Order 4 is also implemented in the butcher tableau, but the rungekutta_4.m implementation + % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it. + if (order == 4) + obj.scheme = @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.scheme = @(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 obj = step(obj) + obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- a/+time/Rungekutta4.m Thu May 02 13:25:14 2019 -0700 +++ b/+time/Rungekutta4.m Wed Aug 07 15:23:42 2019 +0200 @@ -39,7 +39,7 @@ end function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); + obj.v = time.rk.rungekutta_4(obj.v, obj.t, obj.k, obj.F); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- a/+time/Rungekutta4SecondOrder.m Thu May 02 13:25:14 2019 -0700 +++ b/+time/Rungekutta4SecondOrder.m Wed Aug 07 15:23:42 2019 +0200 @@ -99,7 +99,7 @@ end function obj = step(obj) - obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); + obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- a/+time/Rungekutta4proper.m Thu May 02 13:25:14 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,42 +0,0 @@ -classdef Rungekutta4proper < time.Timestepper - properties - F - k - t - v - m - n - end - - - methods - % Timesteps v_t = F(v,t), using RK4 fromt t = t0 with timestep k and initial conditions v = v0 - function obj = Rungekutta4proper(F, k, t0, v0) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.m = length(v0); - obj.n = 0; - end - - function [v,t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function obj = step(obj) - obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - end - end - - - methods (Static) - function k = getTimeStep(lambda) - k = rk4.get_rk4_time_step(lambda); - end - end - -end \ No newline at end of file