Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 591:39554f2de783 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Utux2D scheme
author Martin Almquist <martin.almquist@it.uu.se>
date Mon, 11 Sep 2017 14:12:54 +0200
parents
children 2a2f34778ded
comparison
equal deleted inserted replaced
573:efe2dbf9796e 591:39554f2de783
1 classdef Utux2D < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 grid % Grid
6 order % Order accuracy for the approximation
7 v0 % Initial data
8
9 a % Wave speed a = [a1, a2];
10
11 H % Discrete norm
12 Hi, Hx, Hy, Hxi, Hyi
13
14 % Derivatives
15 Dx, Dy
16
17 % Boundary operators
18 e_w, e_e, e_s, e_n
19
20 D % Total discrete operator
21
22 end
23
24
25 methods
26 function obj = Utux2D(g ,order, opSet, a)
27
28 default_arg('a',1/sqrt(2)*[1, 1]);
29 default_arg('opSet',@sbp.D2Standard);
30 assert(isa(g, 'grid.Cartesian'))
31
32 m = g.size();
33 m_x = m(1);
34 m_y = m(2);
35 m_tot = g.N();
36
37 xlim = g.x{1};
38 ylim = g.x{2};
39 obj.grid = g;
40
41 % Operator sets
42 ops_x = opSet(m_x, xlim, order);
43 ops_y = opSet(m_y, ylim, order);
44 Ix = speye(m_x);
45 Iy = speye(m_y);
46
47 % Norms
48 Hx = ops_x.H;
49 Hy = ops_y.H;
50 Hxi = ops_x.HI;
51 Hyi = ops_y.HI;
52 obj.H = kron(Hx,Hy);
53 obj.Hi = kron(Hxi,Hyi);
54 obj.Hx = kron(Hx,Iy);
55 obj.Hy = kron(Ix,Hy);
56 obj.Hxi = kron(Hxi,Iy);
57 obj.Hyi = kron(Ix,Hyi);
58
59 % Derivatives
60 Dx = ops_x.D1;
61 Dy = ops_y.D1;
62 obj.Dx = kron(Dx,Iy);
63 obj.Dy = kron(Ix,Dy);
64
65 % Boundary operators
66 obj.e_w = kr(ops_x.e_l, Iy);
67 obj.e_e = kr(ops_x.e_r, Iy);
68 obj.e_s = kr(Ix, ops_y.e_l);
69 obj.e_n = kr(Ix, ops_y.e_r);
70
71 obj.m = m;
72 obj.h = [ops_x.h ops_y.h];
73 obj.order = order;
74
75 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
76
77 end
78 % Closure functions return the opertors applied to the own domain to close the boundary
79 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
80 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
81 % type is a string specifying the type of boundary condition if there are several.
82 % data is a function returning the data that should be applied at the boundary.
83 % neighbour_scheme is an instance of Scheme that should be interfaced to.
84 % neighbour_boundary is a string specifying which boundary to interface to.
85 function [closure, penalty] = boundary_condition(obj,boundary,type)
86 default_arg('type','dirichlet');
87
88 sigma = -1; % Scalar penalty parameter
89 switch boundary
90 case {'w','W','west','West'}
91 tau = sigma*obj.a(1)*obj.e_w*obj.Hy;
92 closure = obj.Hi*tau*obj.e_w';
93
94 case {'s','S','south','South'}
95 tau = sigma*pbj.a(2)*obj.e_s*obj.Hx;
96 closure = obj.Hi*tau*obj.e_s';
97 end
98 penalty = -obj.Hi*tau;
99
100 end
101
102 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
103
104 % Get neighbour boundary operator
105 switch neighbour_boundary
106 case {'e','E','east','East'}
107 e_neighbour = neighbour_scheme.e_e;
108 case {'w','W','west','West'}
109 e_neighbour = neighbour_scheme.e_w;
110 case {'n','N','north','North'}
111 e_neighbour = neighbour_scheme.e_n;
112 case {'s','S','south','South'}
113 e_neighbour = neighbour_scheme.e_s;
114 end
115
116 % Upwind coupling
117 sigma_ds = -1; %"Downstream" penalty
118 sigma_us = 0; %"Upstream" penalty
119
120 switch boundary
121 case {'w','W','west','West'}
122 tau = sigma_ds*obj.a(1)*obj.e_w*obj.Hy;
123 closure = obj.Hi*tau*obj.e_w';
124 case {'e','E','east','East'}
125 tau = sigma_us*obj.a(1)*obj.e_e*obj.Hy;
126 closure = obj.Hi*tau*obj.e_e';
127 case {'s','S','south','South'}
128 tau = sigma_ds*obj.a(2)*obj.e_s*obj.Hx;
129 closure = obj.Hi*tau*obj.e_s';
130 case {'n','N','north','North'}
131 tau = sigma_us*obj.a(2)*obj.e_n*obj.Hx;
132 closure = obj.Hi*tau*obj.e_n';
133 end
134 penalty = -obj.Hi*tau*e_neighbour';
135
136 end
137
138 function N = size(obj)
139 N = obj.m;
140 end
141
142 end
143
144 methods(Static)
145 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
146 % and bound_v of scheme schm_v.
147 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
148 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
149 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
150 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
151 end
152 end
153 end