Mercurial > repos > public > sbplib
changeset 1012:1e437c9e5132 feature/advectionRV
Create residual viscosity package +rv and generalize the ResidualViscosity class
- Generalize residual viscosity, by passing user-defined flux and calculating the time derivative outside of the update.
- Create separate RungekuttaRV specifically using interior RV updates
- Separate the artifical dissipation operator from the scheme AdvectionRV1D so that the same scheme can be reused for creating the diff op used by the ResidualViscosity class
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 05 Dec 2018 13:44:10 +0100 |
| parents | e0560bc4fb7d |
| children | eb441fbdf379 |
| files | +rv/+time/RungekuttaInteriorRV.m +rv/+time/rungekuttaRV.m +rv/ResidualViscosity.m +scheme/AdvectionRV1D.m +time/+rk/rungekuttaRV.m +time/RungekuttaRV.m |
| diffstat | 6 files changed, 143 insertions(+), 90 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/RungekuttaInteriorRV.m Wed Dec 05 13:44:10 2018 +0100 @@ -0,0 +1,42 @@ +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 + end + + methods + function obj = RungekuttaInteriorRV(F, k, t0, v0, RV, order) + obj.F = F; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.n = 0; + obj.RV = RV; + % 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); + end + + function [v, t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function state = getState(obj) + [residual, u_t, grad_f] = obj.RV.getResidual(); + state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t); + end + + function obj = step(obj) + obj.v = rv.time.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+rv/+time/rungekuttaRV.m Wed Dec 05 13:44:10 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, 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.getViscosity()); + + % 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.getViscosity()); + 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 Wed Dec 05 13:44:10 2018 +0100 @@ -0,0 +1,65 @@ +% class describing the viscosity +classdef ResidualViscosity < handle + properties + D % 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 controling magnitude of upwind dissipation + Cres % Constant controling magnitude residual dissipation + h % Length scale used for scaling the viscosity. + viscosity % Stores the computed viscosity. + + % Convenice (for verification and plotting) TBD: Decide on if it should be kept. + u_t % Stores the latest approximated time derivative of the solution. + grad_f % Stores the latest approximated gradient of the flux + residual % Stores the computed residual + end + + methods + % TODO: - Consider passing residual normalization as a function handle. + % or choosing a type of normalization on construction. + % Could for example be 1, norm((v-mean(v),inf) or normInfNeighborhood(v) + % but working + % - 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(D, waveSpeed, Cmax, Cres, h, N) + obj.D = D; + obj.waveSpeed = waveSpeed; + obj.h = h; + obj.Cmax = Cmax; + obj.Cres = Cres; + obj.viscosity = zeros(N,1); + obj.u_t = zeros(N,1); + obj.grad_f = zeros(N,1); + obj.residual = zeros(N,1); + end + + function obj = update(obj, v, dvdt) + obj.u_t = dvdt; + obj.grad_f = obj.D(v); + obj.residual = obj.u_t + obj.grad_f; + obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf)); + end + + function [residual, u_t, grad_f] = getResidual(obj) + residual = obj.residual; + u_t = obj.u_t; + grad_f = obj.grad_f; + end + + function viscosity = getViscosity(obj) + viscosity = obj.viscosity; + end + end + % Remove or fix. Should be able to handle values close to zero. Should work in 2d and 3d. + methods (Static) + function R_norm = normInfNeighborhood(v) + n = length(v); + R_norm = zeros(n,1); + R_norm(1,1) = norm(v(1:3), inf); + R_norm(n,1) = norm(v(n-3:n), inf); + for i = 2:n-1 + R_norm(i,1) = norm(v(i-1:i+1), inf); + end + end + end +end \ No newline at end of file
--- a/+scheme/AdvectionRV1D.m Thu Nov 15 13:49:11 2018 -0800 +++ b/+scheme/AdvectionRV1D.m Wed Dec 05 13:44:10 2018 +0100 @@ -8,6 +8,8 @@ Hi % Norm inverse e_l e_r + + D2_visc % Artificial viscosity operator end methods @@ -19,13 +21,14 @@ 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'; - D2 = @(viscosity) ops.Dm*spdiag(viscosity)*ops.Dp-ops.HI*(B*spdiag(viscosity)*ops.Dp); + 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 - % max(abs()) or just abs()? - obj.D = @(viscosity) (-D1 + D2(viscosity) + DissipationOp); + + obj.D = -D1 + DissipationOp; obj.grid = grid; obj.order = order;
--- a/+time/+rk/rungekuttaRV.m Thu Nov 15 13:49:11 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -% Takes one time step of size dt using the rungekutta method -% starting from v_0 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, 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.getViscosity()); - - % 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(u,v,coeffs.c(i)*dt); - k(:,i) = F(u,t+coeffs.c(i)*dt, RV.getViscosity()); - 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(u,v,dt); - v = u; -end
--- a/+time/RungekuttaRV.m Thu Nov 15 13:49:11 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -classdef RungekuttaRV < time.Timestepper - properties - F % RHS of the ODE - k % Time step - t % Time point - v % Solution vector - n % Time level - RV % Residual Viscosity - coeffs % The coefficents used for the RK time integration - end - - methods - function obj = RungekuttaRV(F, k, t0, v0, RV, order) - obj.F = F; - obj.k = k; - obj.t = t0; - obj.v = v0; - obj.n = 0; - obj.RV = RV; - % 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); - end - - function [v, t] = getV(obj) - v = obj.v; - t = obj.t; - end - - function state = getState(obj) - [residual, u_t, grad_f] = obj.RV.getResidual(); - state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t); - end - - function obj = step(obj) - obj.v = time.rk.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs); - obj.t = obj.t + obj.k; - obj.n = obj.n + 1; - % TBD: Add option for updating the residual inside or outside? Decide on best way to do it? - % v_prev = obj.v; - % F = @(v,t)obj.F(v,t,obj.RV.getViscosity()); - % obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F, obj.coeffs); - % obj.RV.update(obj.v,v_prev,obj.k); - % 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