Mercurial > repos > public > sbplib
changeset 1152:010bb2677230 feature/rv
Clean up in +rv/+time. Make the time stepping more efficient by not storing unnessecary properties in the RK-RV time steppers
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Tue, 05 Mar 2019 10:53:34 +0100 |
| parents | 03ecf18d035f |
| children | 635386c073b9 0c906f7ab8bf |
| files | +rv/+time/BdfDerivative.m +rv/+time/RungekuttaExteriorRv.m +rv/+time/RungekuttaExteriorRvBdf.m +rv/+time/RungekuttaInteriorRv.m +rv/+time/rungekuttaRV.m +rv/ResidualViscosity.m |
| diffstat | 6 files changed, 43 insertions(+), 83 deletions(-) [+] |
line wrap: on
line diff
--- a/+rv/+time/BdfDerivative.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/+time/BdfDerivative.m Tue Mar 05 10:53:34 2019 +0100 @@ -2,16 +2,10 @@ 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. + methods 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;
--- a/+rv/+time/RungekuttaExteriorRv.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/+time/RungekuttaExteriorRv.m Tue Mar 05 10:53:34 2019 +0100 @@ -6,22 +6,8 @@ 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 - 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 - viscosity % Total viscosity - residualViscosity % Residual viscosity - firstOrderViscosity % first order viscosity - dvdt % Evaluated time derivative in residual - Df % Evaluated flux in residual end methods @@ -38,8 +24,6 @@ obj.RV = RV; obj.DvDt = DvDt; - obj.dvdt = obj.DvDt(obj.v); - [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = RV.evaluate(obj.v,obj.dvdt); end function [v, t] = getV(obj) @@ -48,15 +32,16 @@ end function state = getState(obj) - state = struct('v', obj.v, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 'residualViscosity', obj.residualViscosity, 'firstOrderViscosity', obj.firstOrderViscosity, 't', obj.t); + 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) - obj.dvdt = obj.DvDt(obj.v); - [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt); - % Fix the viscosity of the RHS function F - F_visc = @(v,t) obj.F(v,t,obj.viscosity); + F_visc = @(v,t) obj.F(v,t,obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v))); 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;
--- a/+rv/+time/RungekuttaExteriorRvBdf.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/+time/RungekuttaExteriorRvBdf.m Tue Mar 05 10:53:34 2019 +0100 @@ -15,13 +15,6 @@ % 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 - viscosity % Total viscosity - residualViscosity % Residual viscosity - firstOrderViscosity % first order viscosity - dvdt % Evaluated time derivative in residual - Df % Evaluated flux in residual end methods function obj = RungekuttaExteriorRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders) @@ -36,22 +29,11 @@ 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.firstOrderViscosity = zeros(size(v0)); - obj.residualViscosity = zeros(size(v0)); - obj.dvdt = zeros(size(v0)); - obj.Df = zeros(size(v0)); end function [v, t] = getV(obj) @@ -60,31 +42,42 @@ end function state = getState(obj) - state = struct('v', obj.v, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 'residualViscosity', obj.residualViscosity, 'firstOrderViscosity', obj.firstOrderViscosity, 't', obj.t); + 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) - % Store current time level and update v_prev - numStoredStages = size(obj.v_prev,2); - if (numStoredStages < obj.upperBdfOrder) + 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]; - numStoredStages = numStoredStages+1; 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 - F_visc = @(v,t) obj.F(v,t,obj.viscosity); + F_visc = @(v,t) obj.F(v, t, 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.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt); - end end end end \ No newline at end of file
--- a/+rv/+time/RungekuttaInteriorRv.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/+time/RungekuttaInteriorRv.m Tue Mar 05 10:53:34 2019 +0100 @@ -8,13 +8,6 @@ coeffs % The coefficents used for the RK time integration RV % Residual Viscosity DvDt % Function for computing the time deriative used for the RV evaluation - - % Convenience properties. Only for plotting - viscosity % Total viscosity - residualViscosity % Residual viscosity - firstOrderViscosity % first order viscosity - dvdt % Evaluated time derivative in residual - Df % Evaluated flux in residual end methods @@ -30,8 +23,6 @@ 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.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt); end function [v, t] = getV(obj) @@ -40,15 +31,17 @@ end function state = getState(obj) - state = struct('v', obj.v, 'dvdt', obj.dvdt, 'Df', obj.Df, 'viscosity', obj.viscosity, 'residualViscosity', obj.residualViscosity, 'firstOrderViscosity', obj.firstOrderViscosity, 't', obj.t); + 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; - obj.dvdt = obj.DvDt(obj.v); - [obj.viscosity, obj.Df, obj.firstOrderViscosity, obj.residualViscosity] = obj.RV.evaluate(obj.v,obj.dvdt); end end end \ No newline at end of file
--- a/+rv/+time/rungekuttaRV.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/+time/rungekuttaRV.m Tue Mar 05 10:53:34 2019 +0100 @@ -7,7 +7,7 @@ % 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))); + k(:,1) = F(v,t,RV.evaluateViscosity(v,DvDt(v))); % Compute the intermediate stages k for i = 2:coeffs.s @@ -15,8 +15,7 @@ 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))); + k(:,i) = F(u,t+coeffs.c(i)*dt, RV.evaluateViscosity(u,DvDt(u))); end % Compute the updated solution as a linear combination @@ -25,6 +24,5 @@ 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
--- a/+rv/ResidualViscosity.m Mon Feb 18 09:00:00 2019 +0100 +++ b/+rv/ResidualViscosity.m Tue Mar 05 10:53:34 2019 +0100 @@ -24,6 +24,12 @@ obj.boundaryIndices = obj.getBoundaryIndices(g); end + function viscosity = evaluateViscosity(obj, v, dvdt) + viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(dvdt + obj.Df(v))./obj.normalization(v)); + % TODO: If the boundary conditions are implemented correctly, this might not be needed. + viscosity(obj.boundaryIndices) = 0; + end + function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt) Df = obj.Df(v); firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v)); @@ -47,15 +53,6 @@ end end - % function f = minmaxDiffNeighborhood(g) - % switch g.D() - % case 1 - % f = @(u)minmaxDiffNeighborhood1d(u); - % case 2 - % f = @(u)minmaxDiffNeighborhood2d(g,u); - % end - % end - function minmaxDiff = minmaxDiffNeighborhood1d(u) umax = movmax(u,3); umin = movmin(u,3);