Mercurial > repos > public > sbplib
annotate +scheme/Utux2D.m @ 605:0f9d20dbb7ce feature/utux2D
Add centered interface coupling in addition to upwind
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 05 Oct 2017 20:21:20 -0700 |
| parents | 2a2f34778ded |
| children | b7b3c11fab4d |
| rev | line source |
|---|---|
| 591 | 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 | |
| 595 | 12 H_x, H_y % Norms in the x and y directions |
| 13 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms | |
| 591 | 14 |
| 15 % Derivatives | |
| 16 Dx, Dy | |
| 17 | |
| 18 % Boundary operators | |
| 19 e_w, e_e, e_s, e_n | |
| 20 | |
| 21 D % Total discrete operator | |
|
605
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
22 |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
23 % String, type of interface coupling |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
24 % Default: 'upwind' |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
25 % Other: 'centered' |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
26 coupling_type |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
27 |
| 591 | 28 |
| 29 end | |
| 30 | |
| 31 | |
| 32 methods | |
|
605
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
33 function obj = Utux2D(g ,order, opSet, a, coupling_type) |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
34 |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
35 default_arg('coupling_type','upwind'); |
| 591 | 36 default_arg('a',1/sqrt(2)*[1, 1]); |
| 37 default_arg('opSet',@sbp.D2Standard); | |
| 38 assert(isa(g, 'grid.Cartesian')) | |
| 39 | |
| 40 m = g.size(); | |
| 41 m_x = m(1); | |
| 42 m_y = m(2); | |
| 43 m_tot = g.N(); | |
| 44 | |
| 595 | 45 xlim = {g.x{1}(1), g.x{1}(end)}; |
| 46 ylim = {g.x{2}(1), g.x{2}(end)}; | |
| 591 | 47 obj.grid = g; |
| 48 | |
| 49 % Operator sets | |
| 50 ops_x = opSet(m_x, xlim, order); | |
| 51 ops_y = opSet(m_y, ylim, order); | |
| 52 Ix = speye(m_x); | |
| 53 Iy = speye(m_y); | |
| 54 | |
| 55 % Norms | |
| 56 Hx = ops_x.H; | |
| 57 Hy = ops_y.H; | |
| 58 Hxi = ops_x.HI; | |
| 59 Hyi = ops_y.HI; | |
| 595 | 60 |
| 61 obj.H_x = Hx; | |
| 62 obj.H_y = Hy; | |
| 591 | 63 obj.H = kron(Hx,Hy); |
| 64 obj.Hi = kron(Hxi,Hyi); | |
| 65 obj.Hx = kron(Hx,Iy); | |
| 66 obj.Hy = kron(Ix,Hy); | |
| 67 obj.Hxi = kron(Hxi,Iy); | |
| 68 obj.Hyi = kron(Ix,Hyi); | |
| 69 | |
| 70 % Derivatives | |
| 71 Dx = ops_x.D1; | |
| 72 Dy = ops_y.D1; | |
| 73 obj.Dx = kron(Dx,Iy); | |
| 74 obj.Dy = kron(Ix,Dy); | |
| 75 | |
| 76 % Boundary operators | |
| 77 obj.e_w = kr(ops_x.e_l, Iy); | |
| 78 obj.e_e = kr(ops_x.e_r, Iy); | |
| 79 obj.e_s = kr(Ix, ops_y.e_l); | |
| 80 obj.e_n = kr(Ix, ops_y.e_r); | |
| 81 | |
| 82 obj.m = m; | |
| 83 obj.h = [ops_x.h ops_y.h]; | |
| 84 obj.order = order; | |
| 595 | 85 obj.a = a; |
|
605
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
86 obj.coupling_type = coupling_type; |
| 591 | 87 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); |
| 88 | |
| 89 end | |
| 90 % Closure functions return the opertors applied to the own domain to close the boundary | |
| 91 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 92 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 93 % type is a string specifying the type of boundary condition if there are several. | |
| 94 % data is a function returning the data that should be applied at the boundary. | |
| 95 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 96 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 97 function [closure, penalty] = boundary_condition(obj,boundary,type) | |
| 98 default_arg('type','dirichlet'); | |
| 99 | |
| 100 sigma = -1; % Scalar penalty parameter | |
| 101 switch boundary | |
| 102 case {'w','W','west','West'} | |
| 595 | 103 tau = sigma*obj.a(1)*obj.e_w*obj.H_y; |
| 591 | 104 closure = obj.Hi*tau*obj.e_w'; |
| 105 | |
| 106 case {'s','S','south','South'} | |
| 595 | 107 tau = sigma*obj.a(2)*obj.e_s*obj.H_x; |
| 591 | 108 closure = obj.Hi*tau*obj.e_s'; |
| 109 end | |
| 110 penalty = -obj.Hi*tau; | |
| 111 | |
| 112 end | |
| 113 | |
| 114 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 115 | |
| 116 % Get neighbour boundary operator | |
| 117 switch neighbour_boundary | |
| 118 case {'e','E','east','East'} | |
| 119 e_neighbour = neighbour_scheme.e_e; | |
| 120 case {'w','W','west','West'} | |
| 121 e_neighbour = neighbour_scheme.e_w; | |
| 122 case {'n','N','north','North'} | |
| 123 e_neighbour = neighbour_scheme.e_n; | |
| 124 case {'s','S','south','South'} | |
| 125 e_neighbour = neighbour_scheme.e_s; | |
| 126 end | |
| 127 | |
|
605
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
128 switch obj.coupling_type |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
129 |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
130 % Upwind coupling (energy dissipation) |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
131 case 'upwind' |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
132 sigma_ds = -1; %"Downstream" penalty |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
133 sigma_us = 0; %"Upstream" penalty |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
134 |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
135 % Energy-preserving coupling (no energy dissipation) |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
136 case 'centered' |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
137 sigma_ds = -1/2; %"Downstream" penalty |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
138 sigma_us = 1/2; %"Upstream" penalty |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
139 |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
140 otherwise |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
141 error(['Interface coupling type ' coupling_type ' is not available.']) |
|
0f9d20dbb7ce
Add centered interface coupling in addition to upwind
Martin Almquist <malmquist@stanford.edu>
parents:
595
diff
changeset
|
142 end |
| 591 | 143 |
| 144 switch boundary | |
| 145 case {'w','W','west','West'} | |
| 595 | 146 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; |
| 591 | 147 closure = obj.Hi*tau*obj.e_w'; |
| 148 case {'e','E','east','East'} | |
| 595 | 149 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; |
| 591 | 150 closure = obj.Hi*tau*obj.e_e'; |
| 151 case {'s','S','south','South'} | |
| 595 | 152 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; |
| 591 | 153 closure = obj.Hi*tau*obj.e_s'; |
| 154 case {'n','N','north','North'} | |
| 595 | 155 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; |
| 591 | 156 closure = obj.Hi*tau*obj.e_n'; |
| 157 end | |
| 158 penalty = -obj.Hi*tau*e_neighbour'; | |
| 159 | |
| 160 end | |
| 161 | |
| 162 function N = size(obj) | |
| 163 N = obj.m; | |
| 164 end | |
| 165 | |
| 166 end | |
| 167 | |
| 168 methods(Static) | |
| 169 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u | |
| 170 % and bound_v of scheme schm_v. | |
| 171 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') | |
| 172 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 173 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 174 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 175 end | |
| 176 end | |
| 177 end |