Mercurial > repos > public > sbplib
changeset 594:a2ddaccf5fd1 feature/utux2D
Merge with better_multiblock_defs
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 18 Sep 2017 20:06:25 +0200 |
| parents | 4422c4476650 (diff) 37948bfe9d79 (current diff) |
| children | 2a2f34778ded |
| files | +multiblock/Def.m |
| diffstat | 2 files changed, 188 insertions(+), 26 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Utux.m Mon Sep 18 19:41:12 2017 +0200 +++ b/+scheme/Utux.m Mon Sep 18 20:06:25 2017 +0200 @@ -2,7 +2,7 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing - x % Grid + grid % Grid order % Order accuracy for the approximation H % Discrete norm @@ -17,41 +17,40 @@ methods - function obj = Utux(m,xlim,order,operator) - default_arg('a',1); - - %Old operators - % [x, h] = util.get_grid(xlim{:},m); - %ops = sbp.Ordinary(m,h,order); - - + function obj = Utux(g ,order, operator) + default_arg('operator','Standard'); + + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + switch operator - case 'NonEquidistant' - ops = sbp.D1Nonequidistant(m,xlim,order); - obj.D1 = ops.D1; +% case 'NonEquidistant' +% ops = sbp.D1Nonequidistant(m,xlim,order); +% obj.D1 = ops.D1; case 'Standard' ops = sbp.D2Standard(m,xlim,order); obj.D1 = ops.D1; - case 'Upwind' - ops = sbp.D1Upwind(m,xlim,order); - obj.D1 = ops.Dm; +% case 'Upwind' +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; otherwise error('Unvalid operator') end - obj.x=ops.x; + + obj.grid = g; - obj.H = ops.H; obj.Hi = ops.HI; obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.D=obj.D1; + obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -61,17 +60,27 @@ % 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','neumann'); - default_arg('data',0); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; + penalty = -obj.Hi*tau; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + end function N = size(obj) @@ -81,9 +90,9 @@ end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Utux2D.m Mon Sep 18 20:06:25 2017 +0200 @@ -0,0 +1,153 @@ +classdef Utux2D < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + grid % Grid + order % Order accuracy for the approximation + v0 % Initial data + + a % Wave speed a = [a1, a2]; + + H % Discrete norm + Hi, Hx, Hy, Hxi, Hyi + + % Derivatives + Dx, Dy + + % Boundary operators + e_w, e_e, e_s, e_n + + D % Total discrete operator + + end + + + methods + function obj = Utux2D(g ,order, opSet, a) + + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + assert(isa(g, 'grid.Cartesian')) + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = g.x{1}; + ylim = g.x{2}; + 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 = kron(Hx,Hy); + obj.Hi = kron(Hxi,Hyi); + obj.Hx = kron(Hx,Iy); + obj.Hy = kron(Ix,Hy); + obj.Hxi = kron(Hxi,Iy); + obj.Hyi = kron(Ix,Hyi); + + % Derivatives + Dx = ops_x.D1; + Dy = ops_y.D1; + obj.Dx = kron(Dx,Iy); + obj.Dy = kron(Ix,Dy); + + % 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); + + obj.m = m; + obj.h = [ops_x.h ops_y.h]; + obj.order = order; + + obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); + + end + % Closure functions return the opertors applied to the own domain 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 doamin. + % 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) + default_arg('type','dirichlet'); + + sigma = -1; % Scalar penalty parameter + switch boundary + case {'w','W','west','West'} + tau = sigma*obj.a(1)*obj.e_w*obj.Hy; + closure = obj.Hi*tau*obj.e_w'; + + case {'s','S','south','South'} + tau = sigma*pbj.a(2)*obj.e_s*obj.Hx; + closure = obj.Hi*tau*obj.e_s'; + end + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + end + + % Upwind coupling + sigma_ds = -1; %"Downstream" penalty + sigma_us = 0; %"Upstream" penalty + + switch boundary + case {'w','W','west','West'} + tau = sigma_ds*obj.a(1)*obj.e_w*obj.Hy; + closure = obj.Hi*tau*obj.e_w'; + case {'e','E','east','East'} + tau = sigma_us*obj.a(1)*obj.e_e*obj.Hy; + closure = obj.Hi*tau*obj.e_e'; + case {'s','S','south','South'} + tau = sigma_ds*obj.a(2)*obj.e_s*obj.Hx; + closure = obj.Hi*tau*obj.e_s'; + case {'n','N','north','North'} + tau = sigma_us*obj.a(2)*obj.e_n*obj.Hx; + closure = obj.Hi*tau*obj.e_n'; + end + penalty = -obj.Hi*tau*e_neighbour'; + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file