Mercurial > repos > public > sbplib
changeset 1034:2676ad79f994 feature/burgers1d
Rename Burgers1D to Burgers1d
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 17 Jan 2019 11:03:22 +0100 |
| parents | 037f203b9bf5 |
| children | 2b9bdb22baec |
| files | +scheme/Burgers1D.m +scheme/Burgers1d |
| diffstat | 2 files changed, 148 insertions(+), 148 deletions(-) [+] |
line wrap: on
line diff
diff -r 037f203b9bf5 -r 2676ad79f994 +scheme/Burgers1D.m --- a/+scheme/Burgers1D.m Thu Jan 17 10:44:12 2019 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,148 +0,0 @@ -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
diff -r 037f203b9bf5 -r 2676ad79f994 +scheme/Burgers1d --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Burgers1d Thu Jan 17 11:03:22 2019 +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