Mercurial > repos > public > sbplib
changeset 1018:eca4d9df9365 feature/advectionRV
Merge with default for latest changes
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Tue, 11 Dec 2018 17:31:28 +0100 |
| parents | 2d7c1333bd6c (diff) 92c3c170e90b (current diff) |
| children | f029b97dbc72 |
| files | +sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m +sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m +sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m +sbp/InterpAWW.m +sbp/InterpMC.m +scheme/Wave.m |
| diffstat | 24 files changed, 1194 insertions(+), 165 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/BdfDerivative.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,31 @@ +classdef BdfDerivative < handle + properties + coefficients + end + methods %TBD: Decide if the BDF order should be passed in construction + % and only a row of coefficients stored based on the order. + % There would still be an implicit dependancy on the number + % of vectors in v_prev and elements in coefficients. + % In addition, dt could be stored, but this would be + % inflexible if different step sizes are employed. + function obj = BdfDerivative() + obj.coefficients = obj.getBdfCoefficients(); + end + % Add asserts here? + 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(6,7); + c(1,1) = 1; c(1,2) = 1; + c(2,1) = 3/2; c(2,2) = 4/2; c(2,3) = -1/2; + c(3,1) = 11/6; c(3,2) = 18/6; c(3,3) = -9/6; c(3,4) = 2/6; + c(4,1) = 25/12; c(4,2) = 48/12; c(4,3) = -36/12; c(4,4) = 16/12; c(4,5) = -3/12; + c(5,1) = 137/60; c(5,2) = 300/60; c(5,3) = -300/60; c(5,4) = 200/60; c(5,5) = -75/60; c(5,6) = 12/60; + c(6,1) = 147/60; c(6,2) = 360/60; c(6,3) = -450/60; c(6,4) = 400/60; c(6,5) = -225/60; c(6,6) = 72/60; c(6,7) = -10/60; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaExteriorRV.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,101 @@ +classdef RungekuttaExteriorRV < 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 + + % Properties related to the residual viscositys + RV % Residual Viscosity operator + viscosity % Viscosity vector + 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. + + % Convenience properties. Only for plotting + residual + dvdt + Df + end + methods + + % TODO: Decide on how to compute dvdt. + function obj = RungekuttaExteriorRV(F, k, t0, v0, RV, DvDt, rkOrder, bdfOrders) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + % Extract the coefficients for the specified rkOrder + % used for the RK updates from the Butcher tableua. + [s,a,b,c] = time.rk.butcherTableau(rkOrder); + obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); + + obj.RV = RV; + % TBD: Decide on if the initialization of the previous stages used by + % the BDF should be done here, or if it should be checked for each + % step taken. + % If it is moved here, then multiple branching stages can be removed in step() + % but this will effectively result in a plotted simulation starting from n = upperBdfOrder. + % In addition, the properties lowerBdfOrder and upperBdfOrder can be removed. + % obj.lowerBdfOrder = bdfOrders.lowerBdfOrder; + % obj.upperBdfOrder = bdfOrders.upperBdfOrder; + % assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6)); + % obj.v_prev = []; + % obj.DvDt = rv.time.BdfDerivative(); + % obj.viscosity = zeros(size(v0)); + % obj.residual = zeros(size(v0)); + % obj.dvdt = zeros(size(v0)); + % obj.Df = zeros(size(v0)); + + % Using the ODE: + obj.DvDt = DvDt; + obj.dvdt = obj.DvDt(obj.v); + [obj.viscosity, obj.Df] = RV.evaluate(obj.v,obj.dvdt); + obj.residual = obj.dvdt + obj.Df; + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t); + end + + function obj = step(obj) + % Store current time level and update v_prev + % numStoredStages = size(obj.v_prev,2); + % if (numStoredStages < obj.upperBdfOrder) + % obj.v_prev = [obj.v, obj.v_prev]; + % numStoredStages = numStoredStages+1; + % else + % obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1); + % obj.v_prev(:,1) = obj.v; + % end + + obj.dvdt = obj.DvDt(obj.v); + [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); + obj.residual = obj.dvdt + obj.Df; + + % Fix the viscosity of the RHS function F + F_visc = @(v,t) obj.F(v,t,obj.viscosity); + obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + + % %Calculate dvdt and evaluate RV for the new time level + % if ((numStoredStages >= obj.lowerBdfOrder) && (numStoredStages <= obj.upperBdfOrder)) + % obj.dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k); + % [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); + % obj.residual = obj.dvdt + obj.Df; + % end + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaInteriorRV.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,54 @@ +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 + viscosity % Viscosity vector + DvDt % Function for computing the time deriative used for the RV evaluation + + residual + dvdt + Df + 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; + obj.dvdt = obj.DvDt(obj.v); + [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); + obj.residual = obj.dvdt + obj.Df; + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + state = struct('v', obj.v, 'residual', obj.residual, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 't', obj.t); + end + + 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; + obj.dvdt = obj.DvDt(obj.v); + [obj.viscosity, obj.Df] = obj.RV.evaluate(obj.v,obj.dvdt); + obj.residual = obj.dvdt + obj.Df; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/rungekuttaRV.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,30 @@ +% 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. +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.evaluate(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 + %RV.update(0.5*(u+v),(u-v)/(coeffs.c(i)*dt)); % Crank-Nicholson for time discretization + k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluate(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 + %RV.update(0.5*(u+v),(u-v)/dt); % Crank-Nicholson for time discretization + v = u; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/ResidualViscosity.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,30 @@ +classdef ResidualViscosity < handle + properties + 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. + end + + methods + % 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(Df, waveSpeed, Cmax, Cres, h, normalization) + default_arg('normalization',@(v)norm(v-mean(v),inf)); + obj.Df = Df; + obj.waveSpeed = waveSpeed; + obj.h = h; + obj.Cmax = Cmax; + obj.Cres = Cres; + obj.normalization = normalization; + end + + function [viscosity, Df] = evaluate(obj, v, dvdt) + Df = obj.Df(v); + viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + Df)/obj.normalization(v)); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/dissipationOperator.m Tue Dec 11 17:31:28 2018 +0100 @@ -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/AdvectionRV1D.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,91 @@ +classdef AdvectionRV1D < scheme.Scheme + properties + grid % Physical grid + order % Order accuracy for the approximation + + D % Non-stabalized scheme operator + H % Discrete norm + Hi % Norm inverse + e_l + e_r + + D2_visc % Artificial viscosity operator + end + + methods + function obj = AdvectionRV1D(grid, operator_type, order, waveSpeed) + m = grid.size(); + lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. + switch operator_type + case 'upwind+' + ops = sbp.D1Upwind(m, lim, order); + D1 = (ops.Dp + ops.Dm)/2; + B = ops.e_r*ops.e_r' - ops.e_l*ops.e_l'; + obj.D2_visc = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp); + % max(abs()) or just abs()? + DissipationOp = spdiag(abs(waveSpeed))*(ops.Dp-ops.Dm)/2; + otherwise + error('Other operator types not yet supported', operator_type); + end + + obj.D = -D1 + DissipationOp; + + obj.grid = grid; + obj.order = order; + + obj.H = ops.H; + obj.Hi = ops.HI; + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; + 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. + % data is a function returning the data that should be applied at the boundary. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','robin'); + default_arg('data',0); + [e, s] = obj.get_boundary_ops(boundary); + switch type + case {'D', 'dirichlet'} + p = s*obj.Hi*e; + closure = p*e'; + otherwise + error('No such boundary condition: type = %s',type); + end + switch class(data) + case 'double' + penalty = s*p*data; + case 'function_handle' + penalty = @(t) s*p*data(t); + otherwise + error('Wierd data argument!') + end + end + + % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e, s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + s = -1; + case 'r' + e = obj.e_r; + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + 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/AdvectionRV2D.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,223 @@ +classdef AdvectionRV2D < scheme.Scheme + properties + grid % Physical grid + order % Order accuracy for the approximation + + D % Non-stabilized scheme operator + H % Discrete norm + H_inv % Norm inverse + halfnorm_inv % Cell array of inverse halfnorm operators + e_l % Cell array of left boundary operators + e_r % Cell array of right boundary operators + d_l % Cell array of left boundary derivative operators + d_r % Cell array of right boundary derivative operators + waveSpeed + end + + methods + function obj = AdvectionRV2D(g, operator_type, order, dissipation, waveSpeed) + if ~isa(g, 'grid.Cartesian') || g.D() ~= 2 + error('Grid must be 2d cartesian'); + end + + obj.grid = g; + obj.order = order; + + % Create cell array of 1D operators. For example D1_1d{1} = D1_x, D1_1d{2} = D1_y. + [Dp_1d, Dm_1d, H_1d, H_inv_1d, d_l_1d, d_r_1d, e_l_1d, e_r_1d, I, DissipationOp_1d] = ... + obj.assemble1DOperators(g, operator_type, order, dissipation); + + %% 2D-operators + % D1 + D1_1d{1} = (Dp_1d{1} + Dp_1d{1})/2; + D1_1d{2} = (Dp_1d{2} + Dp_1d{2})/2; + D1_2d = obj.extendOperatorTo2D(D1_1d, I); + D1 = D1_2d{1} + D1_2d{2}; + % D2 + + Dp_2d = obj.extendOperatorTo2D(Dp_1d, I); + Dm_2d = obj.extendOperatorTo2D(Dm_1d, I); + D2 = @(viscosity) Dm_2d{1}*spdiag(viscosity)*Dp_2d{1} + Dm_2d{2}*spdiag(viscosity)*Dp_2d{2}; + % m = g.size(); + % ind = grid.funcToMatrix(g, 1:g.N()); + % for i = 1:g.D() + % D2_2d{i} = sparse(zeros(g.N())); + % end + % % x-direction + % for i = 1:m(2) + % p = ind(:,i); + % D2_2d{1}(p,p) = @(viscosity) D2_1d{1}(viscosity(p)); + % end + % % y-direction + % for i = 1:m(1) + % p = ind(i,:); + % D2_2d{2}(p,p) = @(viscosity) D2_1d{2}(viscosity(p)); + % end + % D2 = D2_2d{1} + D2_2d{2}; + + obj.d_l = obj.extendOperatorTo2D(d_l_1d, I); + obj.d_r = obj.extendOperatorTo2D(d_r_1d, I); + obj.e_l = obj.extendOperatorTo2D(e_l_1d, I); + obj.e_r = obj.extendOperatorTo2D(e_r_1d, I); + obj.H = kron(H_1d{1},H_1d{2}); + obj.H_inv = kron(H_inv_1d{1},H_inv_1d{2}); + obj.halfnorm_inv = obj.extendOperatorTo2D(H_inv_1d, I); + obj.waveSpeed = waveSpeed; + + % Dissipation operator + switch dissipation + case 'on' + DissOp_2d = obj.extendOperatorTo2D(DissipationOp_1d, I); + DissOp = DissOp_2d{1} + DissOp_2d{2}; + % max(abs()) or just abs()? + obj.D = @(v, viscosity) (-waveSpeed.*D1 + D2(viscosity) + abs(waveSpeed).*DissOp)*v; + case 'off' + obj.D = @(v, viscosity) (-waveSpeed.*D1 + D2(viscosity))*v; + end + 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. + % data is a function returning the data that should be applied at the boundary. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','robin'); + default_arg('data',0); + [e, d, halfnorm_inv, i_b, s] = obj.get_boundary_ops(boundary); + switch type + case {'D', 'dirichlet'} + p = s*halfnorm_inv*e; + closure = @(v,viscosity) p*(v(i_b)); + case {'N', 'nuemann'} + p = s*halfnorm_inv*e; + closure = @(v,viscosity) p*(viscosity(i_b).*d*v); + case {'R', 'robin'} + p = s*halfnorm_inv*e; + closure = @(v, viscosity) p*(obj.waveSpeed(i_b).*v(i_b) - 2*viscosity(i_b).*d*v); + otherwise + error('No such boundary condition: type = %s',type); + end + switch class(data) + case 'double' + penalty = s*p*data; + case 'function_handle' + penalty = @(t) s*p*data(t); + otherwise + error('Wierd data argument!') + end + end + + % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary. + % The right boundary for each coordinate direction is considered the positive boundary + function [e, d, halfnorm_inv, ind_boundary, s] = get_boundary_ops(obj,boundary) + ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N()); + switch boundary + case 'w' + e = obj.e_l{1}; + d = obj.d_l{1}; + halfnorm_inv = obj.halfnorm_inv{1}; + ind_boundary = ind(1,:); + s = -1; + case 'e' + e = obj.e_r{1}; + d = obj.d_r{1}; + halfnorm_inv = obj.halfnorm_inv{1}; + + ind_boundary = ind(end,:); + s = 1; + case 's' + e = obj.e_l{2}; + d = obj.d_l{2}; + halfnorm_inv = obj.halfnorm_inv{2}; + ind_boundary = ind(:,1); + s = -1; + case 'n' + e = obj.e_r{2}; + d = obj.d_r{2}; + halfnorm_inv = obj.halfnorm_inv{2}; + ind_boundary = ind(:,end); + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + 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 + + methods(Static) + function [Dp, Dm, H, Hi, d_l, d_r, e_l, e_r, I, DissipationOp] = assemble1DOperators(g, operator_type, order, dissipation) + dim = g.D(); + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + DissipationOp = cell(dim,1); + for i = 1:dim + switch operator_type + % case 'narrow' + % ops = sbp.D4Variable(g.m(i), g.lim{i}, order); + % D1{i} = ops.D1; + % D2{i} = ops.D2; + % d_l{i} = ops.d1_l'; + % d_r{i} = ops.d1_r'; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = -1*sbp.dissipationOperator(g.m(i), order, ops.HI); + % end + % case 'upwind-' + % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) ops.Dp*spdiag(viscosity)*ops.Dm; + % d_l{i} = ops.e_l'*ops.Dm; + % d_r{i} = ops.e_r'*ops.Dm; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = (ops.Dp-ops.Dm)/2; + % end + case 'upwind+' + ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + Dp{i} = ops.Dp; + Dm{i} = ops.Dm; + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp; + d_l{i} = ops.e_l'*ops.Dp; + d_r{i} = ops.e_r'*ops.Dp; + if (strcmp(dissipation,'on')) + DissipationOp{i} = (ops.Dp-ops.Dm)/2; + end + % case 'upwind+-' + % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) (ops.Dp*spdiag(viscosity)*ops.Dm + ops.Dm*spdiag(viscosity)*ops.Dp)/2; + % d_l{i} = ops.e_l'*D1; + % d_r{i} = ops.e_r'*D1; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = (ops.Dp-ops.Dm)/2; + % end + otherwise + error('Other operator types not yet supported', operator_type); + end + H{i} = ops.H; + Hi{i} = ops.HI; + e_l{i} = ops.e_l; + e_r{i} = ops.e_r; + I{i} = speye(g.m(i)); + end + end + function op_2d = extendOperatorTo2D(op, I) + op_2d{1} = kr(op{1}, I{2}); + op_2d{2} = kr(I{1}, op{2}); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Burgers1D.m Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,148 @@ +classdef Burgers1D < scheme.Scheme + properties + grid % Physical grid + order % Order accuracy for the approximation + + params + + D % Non-stabalized scheme operator + H % Discrete norm + Hi % Norm inverse + e_l + e_r + d_l + d_r + end + + methods + function obj = Burgers1D(grid, pde_form, operator_type, order, dissipation, params) + assert(grid.D == 1); + assert(grid.size() == length(params.eps)); + m = grid.size(); + lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids. + switch operator_type + case 'narrow' + ops = sbp.D4Variable(m, lim, order); + D1 = ops.D1; + D2 = ops.D2; + if (strcmp(dissipation,'on')) + DissipationOp = -1*sbp.dissipationOperator(m, order, ops.HI); + end + d_l = ops.d1_l'; + d_r = ops.d1_r'; + case 'upwind-' + ops = sbp.D1Upwind(m, lim, order); + D1 = (ops.Dp + ops.Dm)/2; + D2 = @(eps) ops.Dp*spdiag(eps)*ops.Dm; + if (strcmp(dissipation,'on')) + DissipationOp = (ops.Dp-ops.Dm)/2; + end + d_l = ops.e_l'*ops.Dm; + d_r = ops.e_r'*ops.Dm; + case 'upwind+' + ops = sbp.D1Upwind(m, lim, order); + D1 = (ops.Dp + ops.Dm)/2; + D2 = @(eps) ops.Dm*spdiag(eps)*ops.Dp; + if (strcmp(dissipation,'on')) + DissipationOp = (ops.Dp-ops.Dm)/2; + end + d_l = ops.e_l'*ops.Dp; + d_r = ops.e_r'*ops.Dp; + case 'upwind+-' + ops = sbp.D1Upwind(m, lim, order); + D1 = (ops.Dp + ops.Dm)/2; + D2 = @(eps) (ops.Dp*spdiag(eps)*ops.Dm + ops.Dm*spdiag(eps)*ops.Dp)/2; + if (strcmp(dissipation,'on')) + DissipationOp = (ops.Dp-ops.Dm)/2; + end + d_l = ops.e_l'*D1; + d_r = ops.e_r'*D1; + otherwise + error('Other operator types not yet supported', operator_type); + end + + switch pde_form + case 'skew-symmetric' + if (strcmp(dissipation,'on')) + D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; + else + D = @(v, viscosity) - 1/3*D1*v.^2 + (-1/3*v.*D1 + D2(params.eps + viscosity))*v; + end + case 'conservative' + if (strcmp(dissipation,'on')) + D = @(v, viscosity) -1/2*D1*v.^2 + (D2(params.eps + viscosity) + max(abs(v))*DissipationOp)*v; + else + D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; + end + otherwise + error('Not supported', pde_form); + end + + obj.grid = grid; + obj.order = order; + obj.params = params; + + obj.D = D; + obj.H = ops.H; + obj.Hi = ops.HI; + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; + obj.d_l = d_l; + obj.d_r = d_r; + 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. + % data is a function returning the data that should be applied at the boundary. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','robin'); + default_arg('data',0); + [e, d, i_b, s] = obj.get_boundary_ops(boundary); + switch type + % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary + case {'R','robin'} + p = s*obj.Hi*e; + closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3)*(v(i_b)) - ((obj.params.eps(i_b) + viscosity(i_b))*d*v)); + switch class(data) + case 'double' + penalty = s*p*data; + case 'function_handle' + penalty = @(t) s*p*data(t); + otherwise + error('Wierd data argument!') + end + otherwise + error('No such boundary condition: type = %s',type); + end + end + + % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e, d, i_b, s] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d = obj.d_l; + i_b = 1; + s = -1; + case 'r' + e = obj.e_r; + d = obj.d_r; + i_b = length(e); + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + 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 Tue Dec 11 17:31:28 2018 +0100 @@ -0,0 +1,215 @@ +classdef Burgers2D < scheme.Scheme + properties + grid % Physical grid + order % Order accuracy for the approximation + + D % Non-stabilized scheme operator + H % Discrete norm + H_inv % Norm inverse + halfnorm_inv % Cell array halfnorm operators + e_l % Cell array of left boundary operators + e_r % Cell array of right boundary operators + d_l % Cell array of left boundary derivative operators + d_r % Cell array of right boundary derivative operators + end + + methods + function obj = Burgers2D(g, operator_type, order, dissipation) + if ~isa(g, 'grid.Cartesian') || g.D() ~= 2 + error('Grid must be 2d cartesian'); + end + + obj.grid = g; + obj.order = order; + + % Create cell array of 1D operators. For example D1_1d{1} = D1_x, D1_1d{2} = D1_y. + [Dp_1d, Dm_1d, H_1d, H_inv_1d, d_l_1d, d_r_1d, e_l_1d, e_r_1d, I, DissipationOp_1d] = ... + obj.assemble1DOperators(g, operator_type, order, dissipation); + + %% 2D-operators + % D1 + D1_1d{1} = (Dp_1d{1} + Dp_1d{1})/2; + D1_1d{2} = (Dp_1d{2} + Dp_1d{2})/2; + D1_2d = obj.extendOperatorTo2D(D1_1d, I); + D1 = D1_2d{1} + D1_2d{2}; + % D2 + + Dp_2d = obj.extendOperatorTo2D(Dp_1d, I); + Dm_2d = obj.extendOperatorTo2D(Dm_1d, I); + D2 = @(viscosity) Dm_2d{1}*spdiag(viscosity)*Dp_2d{1} + Dm_2d{2}*spdiag(viscosity)*Dp_2d{2}; + % m = g.size(); + % ind = grid.funcToMatrix(g, 1:g.N()); + % for i = 1:g.D() + % D2_2d{i} = sparse(zeros(g.N())); + % end + % % x-direction + % for i = 1:m(2) + % p = ind(:,i); + % D2_2d{1}(p,p) = @(viscosity) D2_1d{1}(viscosity(p)); + % end + % % y-direction + % for i = 1:m(1) + % p = ind(i,:); + % D2_2d{2}(p,p) = @(viscosity) D2_1d{2}(viscosity(p)); + % end + % D2 = D2_2d{1} + D2_2d{2}; + + obj.d_l = obj.extendOperatorTo2D(d_l_1d, I); + obj.d_r = obj.extendOperatorTo2D(d_r_1d, I); + obj.e_l = obj.extendOperatorTo2D(e_l_1d, I); + obj.e_r = obj.extendOperatorTo2D(e_r_1d, I); + obj.H = kron(H_1d{1},H_1d{2}); + obj.H_inv = kron(H_inv_1d{1},H_inv_1d{2}); + obj.halfnorm_inv = obj.extendOperatorTo2D(H_inv_1d, I); + + % Dissipation operator + switch dissipation + case 'on' + DissOp_2d = obj.extendOperatorTo2D(DissipationOp_1d, I); + DissOp = DissOp_2d{1} + DissOp_2d{2}; + obj.D = @(v, viscosity) -1/2*D1*v.^2 + (D2(viscosity) + max(abs(v))*DissOp)*v; + case 'off' + obj.D = @(v, viscosity) -1/2*D1*v.^2 + D2(viscosity)*v; + end + 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. + % data is a function returning the data that should be applied at the boundary. + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','robin'); + default_arg('data',0); + [e, d, halfnorm_inv, i_b, s] = obj.get_boundary_ops(boundary); + switch type + % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary + case {'R','robin'} + p = s*halfnorm_inv*e; + closure = @(v, viscosity) p*(((v(i_b)-s*abs(v(i_b)))/3).*(v(i_b)) - ((viscosity(i_b)).*d*v)); + switch class(data) + case 'double' + penalty = s*p*data; + case 'function_handle' + penalty = @(t) s*p*data(t); + otherwise + error('Wierd data argument!') + end + otherwise + error('No such boundary condition: type = %s',type); + end + end + + % Ruturns the boundary ops, half-norm, boundary indices and sign for the boundary specified by the string boundary. + % The right boundary for each coordinate direction is considered the positive boundary + function [e, d, halfnorm_inv, ind_boundary, s] = get_boundary_ops(obj,boundary) + ind = grid.funcToMatrix(obj.grid, 1:obj.grid.N()); + switch boundary + case 'w' + e = obj.e_l{1}; + d = obj.d_l{1}; + halfnorm_inv = obj.halfnorm_inv{1}; + ind_boundary = ind(1,:); + s = -1; + case 'e' + e = obj.e_r{1}; + d = obj.d_r{1}; + halfnorm_inv = obj.halfnorm_inv{1}; + + ind_boundary = ind(end,:); + s = 1; + case 's' + e = obj.e_l{2}; + d = obj.d_l{2}; + halfnorm_inv = obj.halfnorm_inv{2}; + ind_boundary = ind(:,1); + s = -1; + case 'n' + e = obj.e_r{2}; + d = obj.d_r{2}; + halfnorm_inv = obj.halfnorm_inv{2}; + ind_boundary = ind(:,end); + s = 1; + otherwise + error('No such boundary: boundary = %s',boundary); + 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 + + methods(Static) + function [Dp, Dm, H, Hi, d_l, d_r, e_l, e_r, I, DissipationOp] = assemble1DOperators(g, operator_type, order, dissipation) + dim = g.D(); + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + DissipationOp = cell(dim,1); + for i = 1:dim + switch operator_type + % case 'narrow' + % ops = sbp.D4Variable(g.m(i), g.lim{i}, order); + % D1{i} = ops.D1; + % D2{i} = ops.D2; + % d_l{i} = ops.d1_l'; + % d_r{i} = ops.d1_r'; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = -1*sbp.dissipationOperator(g.m(i), order, ops.HI); + % end + % case 'upwind-' + % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) ops.Dp*spdiag(viscosity)*ops.Dm; + % d_l{i} = ops.e_l'*ops.Dm; + % d_r{i} = ops.e_r'*ops.Dm; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = (ops.Dp-ops.Dm)/2; + % end + case 'upwind+' + ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + Dp{i} = ops.Dp; + Dm{i} = ops.Dm; + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp; + d_l{i} = ops.e_l'*ops.Dp; + d_r{i} = ops.e_r'*ops.Dp; + if (strcmp(dissipation,'on')) + DissipationOp{i} = (ops.Dp-ops.Dm)/2; + end + % case 'upwind+-' + % ops = sbp.D1Upwind(g.m(i), g.lim{i}, order); + % D1{i} = (ops.Dp + ops.Dm)/2; + % D2{i} = @(viscosity) (ops.Dp*spdiag(viscosity)*ops.Dm + ops.Dm*spdiag(viscosity)*ops.Dp)/2; + % d_l{i} = ops.e_l'*D1; + % d_r{i} = ops.e_r'*D1; + % if (strcmp(dissipation,'on')) + % DissipationOp{i} = (ops.Dp-ops.Dm)/2; + % end + otherwise + error('Other operator types not yet supported', operator_type); + end + H{i} = ops.H; + Hi{i} = ops.HI; + e_l{i} = ops.e_l; + e_r{i} = ops.e_r; + I{i} = speye(g.m(i)); + end + end + function op_2d = extendOperatorTo2D(op, I) + op_2d{1} = kr(op{1}, I{2}); + op_2d{2} = kr(I{1}, op{2}); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/butcherTableau.m Tue Dec 11 17:31:28 2018 +0100 @@ -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 Tue Dec 11 17:31:28 2018 +0100 @@ -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 Tue Dec 11 17:31:28 2018 +0100 @@ -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 Tue Dec 11 17:31:28 2018 +0100 @@ -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 Tue Dec 11 17:31:28 2018 +0100 @@ -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 Mon Dec 10 08:29:41 2018 +0000 +++ /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 Mon Dec 10 08:29:41 2018 +0000 +++ /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 Mon Dec 10 08:29:41 2018 +0000 +++ /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 Mon Dec 10 08:29:41 2018 +0000 +++ /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 Mon Dec 10 08:29:41 2018 +0000 +++ b/+time/Rk4SecondOrderNonlin.m Tue Dec 11 17:31:28 2018 +0100 @@ -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 Tue Dec 11 17:31:28 2018 +0100 @@ -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 Mon Dec 10 08:29:41 2018 +0000 +++ b/+time/Rungekutta4.m Tue Dec 11 17:31:28 2018 +0100 @@ -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 Mon Dec 10 08:29:41 2018 +0000 +++ b/+time/Rungekutta4SecondOrder.m Tue Dec 11 17:31:28 2018 +0100 @@ -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 Mon Dec 10 08:29:41 2018 +0000 +++ /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