Mercurial > repos > public > sbplib
changeset 1149:1fe48cbd379a feature/rv
Change Burgers2d to inviscid formulation. Rewrite to use opSets and fix the implementation of the Dirichlet conditions.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 24 Jan 2019 09:05:44 +0100 |
| parents | 0a5503a08a36 |
| children | 922996695952 |
| files | +scheme/Burgers2d.m |
| diffstat | 1 files changed, 113 insertions(+), 169 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Burgers2d.m Thu Jan 24 09:01:10 2019 +0100 +++ b/+scheme/Burgers2d.m Thu Jan 24 09:05:44 2019 +0100 @@ -1,134 +1,147 @@ -classdef Burgers2D < scheme.Scheme +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 + H_x, H_y % Norms in the x and y directions + Hi % Kroneckered norm inverse + % Boundary operators + e_w, e_e, e_s, e_n end methods - function obj = Burgers2D(g, operator_type, order, dissipation) - if ~isa(g, 'grid.Cartesian') || g.D() ~= 2 - error('Grid must be 2d cartesian'); + function obj = Burgers2d(g, order, pde_form, fluxSplitting, opSet) + default_arg('opSet',@sbp.D2Standard); + default_arg('fluxSplitting',{@(v)max(abs(v)),@(v)max(abs(v))}); + assertType(g, 'grid.Cartesian'); + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + obj.grid = g; + + % Operator sets + ops_x = opSet(m_x, xlim, order); + ops_y = opSet(m_y, ylim, order); + Ix = speye(m_x); + Iy = speye(m_y); + + % Norms + Hx = ops_x.H; + Hy = ops_y.H; + Hxi = ops_x.HI; + Hyi = ops_y.HI; + + obj.H_x = Hx; + obj.H_y = Hy; + obj.H = kron(Hx,Hy); + obj.Hi = kron(Hxi,Hyi); + + % Derivatives + if (isequal(opSet,@sbp.D1Upwind)) + Dx = kron((ops_x.Dp + ops_x.Dm)/2,Iy); + Dy = kron(Ix,(ops_y.Dp + ops_y.Dm)/2); + DissOpx = kron((ops_x.Dm - ops_x.Dp)/2,Iy); + DissOpy = kron(Ix,(ops_y.Dm - ops_y.Dp)/2); + D1 = Dx + Dy; + switch pde_form + case 'skew-symmetric' + switch length(fluxSplitting) + case 1 + DissOp = DissOpx + DissOpy; + obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v); + case 2 + obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); + end + case 'conservative' + switch length(fluxSplitting) + case 1 + DissOp = DissOpx + DissOpy; + obj.D = @(v) -(1/2*D1*v.*v + fluxSplitting{1}(v)*DissOp*v); + case 2 + obj.D = @(v) -(1/2*D1*v.*v + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v); + end + + end + else + Dx = kron(ops_x.D1,Iy); + Dy = kron(Ix,ops_y.D1); + D1 = Dx + Dy; + switch pde_form + case 'skew-symmetric' + obj.D = @(v) -(1/3*D1*v.*v + 1/3*spdiag(v)*D1*v); + case 'conservative' + obj.D = @(v) -1/2*D1*v.*v; + end 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 + % Boundary operators + obj.e_w = kr(ops_x.e_l, Iy); + obj.e_e = kr(ops_x.e_r, Iy); + obj.e_s = kr(Ix, ops_y.e_l); + obj.e_n = kr(Ix, ops_y.e_r); - 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 + obj.order = order; 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); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); + [e, H_b, index, 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 + % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3 + % with +- at left/right boundaries in each coordinate direction + case {'D', 'd', 'dirichlet', 'Dirichlet'} + + magnitude = 2/3; + tau = @(v) s*magnitude*obj.Hi*e*H_b*spdiag((v(index)-s*abs(v(index)))/2); + closure = @(v) tau(v)*v(index); + penalty = @(v) -tau(v); + + + % tau = s*e*H_b; + % closure = @(v) (obj.Hi*tau)*(((v(index)-s*abs(v(index)))/3).*v(index)); + % penalty = -obj.Hi*tau; 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) + function [e, H_b, index, 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,:); + case {'w', 'W', 'west', 'West'} + e = obj.e_w; + H_b = obj.H_y; + index = 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,:); + case {'e', 'E', 'east', 'East'} + e = obj.e_e; + H_b = obj.H_y; + index = 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); + case {'s', 'S', 'south', 'South'} + e = obj.e_s; + H_b = obj.H_x; + index = 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); + case {'n', 'N', 'north', 'North'} + e = obj.e_n; + H_b = obj.H_x; + index = ind(:,end); s = 1; otherwise error('No such boundary: boundary = %s',boundary); @@ -143,73 +156,4 @@ 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