Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 605:0f9d20dbb7ce feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
597:8a1c36bbb977 605:0f9d20dbb7ce
17 17
18 % Boundary operators 18 % Boundary operators
19 e_w, e_e, e_s, e_n 19 e_w, e_e, e_s, e_n
20 20
21 D % Total discrete operator 21 D % Total discrete operator
22
23 % String, type of interface coupling
24 % Default: 'upwind'
25 % Other: 'centered'
26 coupling_type
27
22 28
23 end 29 end
24 30
25 31
26 methods 32 methods
27 function obj = Utux2D(g ,order, opSet, a) 33 function obj = Utux2D(g ,order, opSet, a, coupling_type)
28 34
35 default_arg('coupling_type','upwind');
29 default_arg('a',1/sqrt(2)*[1, 1]); 36 default_arg('a',1/sqrt(2)*[1, 1]);
30 default_arg('opSet',@sbp.D2Standard); 37 default_arg('opSet',@sbp.D2Standard);
31 assert(isa(g, 'grid.Cartesian')) 38 assert(isa(g, 'grid.Cartesian'))
32 39
33 m = g.size(); 40 m = g.size();
74 81
75 obj.m = m; 82 obj.m = m;
76 obj.h = [ops_x.h ops_y.h]; 83 obj.h = [ops_x.h ops_y.h];
77 obj.order = order; 84 obj.order = order;
78 obj.a = a; 85 obj.a = a;
86 obj.coupling_type = coupling_type;
79 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); 87 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
80 88
81 end 89 end
82 % Closure functions return the opertors applied to the own domain to close the boundary 90 % Closure functions return the opertors applied to the own domain to close the boundary
83 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 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.
115 e_neighbour = neighbour_scheme.e_n; 123 e_neighbour = neighbour_scheme.e_n;
116 case {'s','S','south','South'} 124 case {'s','S','south','South'}
117 e_neighbour = neighbour_scheme.e_s; 125 e_neighbour = neighbour_scheme.e_s;
118 end 126 end
119 127
120 % Upwind coupling 128 switch obj.coupling_type
121 sigma_ds = -1; %"Downstream" penalty 129
122 sigma_us = 0; %"Upstream" penalty 130 % Upwind coupling (energy dissipation)
131 case 'upwind'
132 sigma_ds = -1; %"Downstream" penalty
133 sigma_us = 0; %"Upstream" penalty
134
135 % Energy-preserving coupling (no energy dissipation)
136 case 'centered'
137 sigma_ds = -1/2; %"Downstream" penalty
138 sigma_us = 1/2; %"Upstream" penalty
139
140 otherwise
141 error(['Interface coupling type ' coupling_type ' is not available.'])
142 end
123 143
124 switch boundary 144 switch boundary
125 case {'w','W','west','West'} 145 case {'w','W','west','West'}
126 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; 146 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y;
127 closure = obj.Hi*tau*obj.e_w'; 147 closure = obj.Hi*tau*obj.e_w';