Mercurial > repos > public > sbplib
view +scheme/Burgers1D.m @ 831:d0934d1143b7 feature/burgers1d
Fix bug in initialization of differential operators
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Tue, 11 Sep 2018 13:24:08 +0200 |
parents | fae41958af4f |
children | 5573913a0949 |
line wrap: on
line source
classdef Burgers1D < scheme.Scheme properties m % Number of points in each direction, possibly a vector h % Grid spacing x % Grid order % Order accuracy for the approximation params D % Non-stabalized scheme operator M % Derivative norm H % Discrete norm Hi % Norm inverse e_l e_r d_l d_r end methods 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(params.eps + viscosity)*v; case 'conservative' D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v; end obj.m = m; obj.h = h; obj.order = order; obj.x = x; obj.params = params; %% 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; obj.e_l = ops.e_l; obj.e_r = ops.e_r; obj.d_l = ops.d1_l; obj.d_r = ops.d1_r; end % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors 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. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type,data) default_arg('type','robin'); default_arg('data',0); [e, s, d] = 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*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*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 % Unknown, boundary condition otherwise error('No such boundary condition: type = %s',type); end end % Ruturns the boundary ops and sign for the boundary specified by the string boundary. % The right boundary is considered the positive boundary function [e, s, d] = get_boundary_ops(obj,boundary) switch boundary case 'l' e = obj.e_l; s = -1; d = obj.d_l; case 'r' e = obj.e_r; s = 1; d = obj.d_r; 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.m; end end end