Mercurial > repos > public > sbplib
diff +scheme/Burgers1D.m @ 815:fae41958af4f feature/burgers1d
Add support for artificial viscosity to the 1d burgers scheme.
- Add support for artificial viscosity by parametrizing the discretization operator and boundary closure on the viscosity.
- Add a time stepper which evaluates and updates the residual viscosity of the solution.
| author | Vidar Stiernstrom <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 06 Sep 2018 12:43:51 +0200 |
| parents | 3a5e635a93fd |
| children | d0934d1143b7 |
line wrap: on
line diff
--- a/+scheme/Burgers1D.m Mon Sep 03 14:50:27 2018 +0200 +++ b/+scheme/Burgers1D.m Thu Sep 06 12:43:51 2018 +0200 @@ -4,6 +4,7 @@ h % Grid spacing x % Grid order % Order accuracy for the approximation + params D % Non-stabalized scheme operator M % Derivative norm @@ -13,21 +14,38 @@ e_r d_l d_r - params % Parameters for the coefficient matrices end methods - function obj = Burgers1D(m, order, xlim, params) - [x, h] = util.get_grid(xlim{:},m); - ops = sbp.D2Variable(m, xlim, order); + function obj = Burgers1D(pde_form, operator_type, order, m, lim, params) + [x, h] = util.get_grid(lim{:},m); + default_arg('pde_form','skew-symmetric'); + default_arg('operator_type','narrow'); + + switch operator_type + case 'narrow' + ops = sbp.D2Variable(m, lim, order); + D1 = ops.D1; + D2 = ops.D2; + otherwise + error('Other operator types not yet supported', operator_type); + end + + switch pde_form + case 'skew-symmetric' + D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(obj.params.eps + viscosity)*v; + case 'conservative' + D = @(v, viscosity) -1/2*D1*v.^2 + D2(obj.params.eps + viscosity)*v; + end obj.m = m; obj.h = h; obj.order = order; obj.x = x; + obj.params = params; - D1 = ops.D1; - obj.D = @(v)(-1/3*v.*D1*v - 1/3*D1*v.^2 + ops.D2(params.eps)*v); + %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity. + obj.D = D; obj.M = ops.M; obj.H = ops.H; obj.Hi = ops.HI; @@ -35,7 +53,6 @@ obj.e_r = ops.e_r; obj.d_l = ops.d1_l; obj.d_r = ops.d1_r; - obj.params = params; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -53,7 +70,7 @@ % 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) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*obj.params.eps*d'*v); + closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v); switch class(data) case 'double' penalty = s*p*data;