Mercurial > repos > public > sbplib
annotate +scheme/Wave2d.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 17 Jan 2019 10:25:06 +0100 |
| parents | 706d1c2b4199 |
| children | 78db023a7fe3 c12b84fe9b00 |
| rev | line source |
|---|---|
|
55
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 classdef Wave2d < scheme.Scheme |
| 0 | 2 properties |
| 3 m % Number of points in each direction, possibly a vector | |
| 4 h % Grid spacing | |
| 5 x,y % Grid | |
| 6 X,Y % Values of x and y for each grid point | |
| 7 order % Order accuracy for the approximation | |
| 8 | |
| 9 D % non-stabalized scheme operator | |
| 10 M % Derivative norm | |
| 11 alpha | |
| 12 | |
| 13 H % Discrete norm | |
| 14 Hi | |
| 15 H_x, H_y % Norms in the x and y directions | |
| 16 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 17 Hi_x, Hi_y | |
| 18 Hix, Hiy | |
| 19 e_w, e_e, e_s, e_n | |
| 20 d1_w, d1_e, d1_s, d1_n | |
| 21 gamm_x, gamm_y | |
| 22 end | |
| 23 | |
| 24 methods | |
|
55
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
25 function obj = Wave2d(m,lim,order,alpha) |
|
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
26 default_arg('alpha',1); |
|
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
27 |
|
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
28 xlim = lim{1}; |
|
a8ed986fcf57
Minor renaming and clean up in 2d wave schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
29 ylim = lim{2}; |
| 0 | 30 |
| 31 if length(m) == 1 | |
| 32 m = [m m]; | |
| 33 end | |
| 34 | |
| 35 m_x = m(1); | |
| 36 m_y = m(2); | |
| 37 | |
| 38 [x, h_x] = util.get_grid(xlim{:},m_x); | |
| 39 [y, h_y] = util.get_grid(ylim{:},m_y); | |
| 40 | |
| 41 ops_x = sbp.Ordinary(m_x,h_x,order); | |
| 42 ops_y = sbp.Ordinary(m_y,h_y,order); | |
| 43 | |
| 44 I_x = speye(m_x); | |
| 45 I_y = speye(m_y); | |
| 46 | |
| 47 D2_x = sparse(ops_x.derivatives.D2); | |
| 48 H_x = sparse(ops_x.norms.H); | |
| 49 Hi_x = sparse(ops_x.norms.HI); | |
| 50 M_x = sparse(ops_x.norms.M); | |
| 51 e_l_x = sparse(ops_x.boundary.e_1); | |
| 52 e_r_x = sparse(ops_x.boundary.e_m); | |
| 53 d1_l_x = sparse(ops_x.boundary.S_1); | |
| 54 d1_r_x = sparse(ops_x.boundary.S_m); | |
| 55 | |
| 56 D2_y = sparse(ops_y.derivatives.D2); | |
| 57 H_y = sparse(ops_y.norms.H); | |
| 58 Hi_y = sparse(ops_y.norms.HI); | |
| 59 M_y = sparse(ops_y.norms.M); | |
| 60 e_l_y = sparse(ops_y.boundary.e_1); | |
| 61 e_r_y = sparse(ops_y.boundary.e_m); | |
| 62 d1_l_y = sparse(ops_y.boundary.S_1); | |
| 63 d1_r_y = sparse(ops_y.boundary.S_m); | |
| 64 | |
| 65 D2 = kr(D2_x, I_y) + kr(I_x, D2_y); | |
| 66 obj.H = kr(H_x,H_y); | |
| 67 obj.Hx = kr(H_x,I_y); | |
| 68 obj.Hy = kr(I_x,H_y); | |
| 69 obj.Hix = kr(Hi_x,I_y); | |
| 70 obj.Hiy = kr(I_x,Hi_y); | |
| 71 obj.Hi = kr(Hi_x,Hi_y); | |
| 72 obj.M = kr(M_x,H_y)+kr(H_x,M_y); | |
| 73 obj.e_w = kr(e_l_x,I_y); | |
| 74 obj.e_e = kr(e_r_x,I_y); | |
| 75 obj.e_s = kr(I_x,e_l_y); | |
| 76 obj.e_n = kr(I_x,e_r_y); | |
| 77 obj.d1_w = kr(d1_l_x,I_y); | |
| 78 obj.d1_e = kr(d1_r_x,I_y); | |
| 79 obj.d1_s = kr(I_x,d1_l_y); | |
| 80 obj.d1_n = kr(I_x,d1_r_y); | |
| 81 | |
| 82 obj.m = m; | |
| 83 obj.h = [h_x h_y]; | |
| 84 obj.order = order; | |
| 85 | |
| 86 obj.alpha = alpha; | |
| 87 obj.D = alpha*D2; | |
| 88 obj.x = x; | |
| 89 obj.y = y; | |
| 90 obj.X = kr(x,ones(m_y,1)); | |
| 91 obj.Y = kr(ones(m_x,1),y); | |
| 92 | |
| 93 obj.gamm_x = h_x*ops_x.borrowing.M.S; | |
| 94 obj.gamm_y = h_y*ops_y.borrowing.M.S; | |
| 95 end | |
| 96 | |
| 97 | |
| 98 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 99 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 100 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 101 % type is a string specifying the type of boundary condition if there are several. | |
| 102 % data is a function returning the data that should be applied at the boundary. | |
| 103 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 104 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 105 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
| 106 default_arg('type','neumann'); | |
| 107 default_arg('data',0); | |
| 108 | |
| 109 [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 110 | |
| 111 switch type | |
| 112 % Dirichlet boundary condition | |
| 113 case {'D','d','dirichlet'} | |
| 114 alpha = obj.alpha; | |
| 115 | |
| 116 % tau1 < -alpha^2/gamma | |
| 117 tuning = 1.1; | |
| 118 tau1 = -tuning*alpha/gamm; | |
| 119 tau2 = s*alpha; | |
| 120 | |
| 121 p = tau1*e + tau2*d; | |
| 122 | |
| 123 closure = halfnorm_inv*p*e'; | |
| 124 | |
| 125 pp = halfnorm_inv*p; | |
| 126 switch class(data) | |
| 127 case 'double' | |
| 128 penalty = pp*data; | |
| 129 case 'function_handle' | |
| 130 penalty = @(t)pp*data(t); | |
| 131 otherwise | |
| 132 error('Wierd data argument!') | |
| 133 end | |
| 134 | |
| 135 | |
| 136 % Neumann boundary condition | |
| 137 case {'N','n','neumann'} | |
| 138 alpha = obj.alpha; | |
| 139 tau1 = -s*alpha; | |
| 140 tau2 = 0; | |
| 141 tau = tau1*e + tau2*d; | |
| 142 | |
| 143 closure = halfnorm_inv*tau*d'; | |
| 144 | |
| 145 pp = halfnorm_inv*tau; | |
| 146 switch class(data) | |
| 147 case 'double' | |
| 148 penalty = pp*data; | |
| 149 case 'function_handle' | |
| 150 penalty = @(t)pp*data(t); | |
| 151 otherwise | |
| 152 error('Wierd data argument!') | |
| 153 end | |
| 154 | |
| 155 % Unknown, boundary condition | |
| 156 otherwise | |
| 157 error('No such boundary condition: type = %s',type); | |
| 158 end | |
| 159 end | |
| 160 | |
|
946
706d1c2b4199
Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents:
910
diff
changeset
|
161 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) |
| 0 | 162 % u denotes the solution in the own domain |
| 163 % v denotes the solution in the neighbour domain | |
| 164 [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 165 [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 166 | |
| 167 tuning = 1.1; | |
| 168 | |
| 169 alpha_u = obj.alpha; | |
| 170 alpha_v = neighbour_scheme.alpha; | |
| 171 | |
| 172 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) | |
| 173 | |
| 174 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; | |
| 175 tau2 = s_u*1/2*alpha_u; | |
| 176 sig1 = s_u*(-1/2); | |
| 177 sig2 = 0; | |
| 178 | |
| 179 tau = tau1*e_u + tau2*d_u; | |
| 180 sig = sig1*e_u + sig2*d_u; | |
| 181 | |
| 182 closure = halfnorm_inv*( tau*e_u' + sig*alpha_u*d_u'); | |
| 183 penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); | |
| 184 end | |
| 185 | |
| 186 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 187 % The right boundary is considered the positive boundary | |
| 188 function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) | |
| 189 switch boundary | |
| 190 case 'w' | |
| 191 e = obj.e_w; | |
| 192 d = obj.d1_w; | |
| 193 s = -1; | |
| 194 gamm = obj.gamm_x; | |
| 195 halfnorm_inv = obj.Hix; | |
| 196 case 'e' | |
| 197 e = obj.e_e; | |
| 198 d = obj.d1_e; | |
| 199 s = 1; | |
| 200 gamm = obj.gamm_x; | |
| 201 halfnorm_inv = obj.Hix; | |
| 202 case 's' | |
| 203 e = obj.e_s; | |
| 204 d = obj.d1_s; | |
| 205 s = -1; | |
| 206 gamm = obj.gamm_y; | |
| 207 halfnorm_inv = obj.Hiy; | |
| 208 case 'n' | |
| 209 e = obj.e_n; | |
| 210 d = obj.d1_n; | |
| 211 s = 1; | |
| 212 gamm = obj.gamm_y; | |
| 213 halfnorm_inv = obj.Hiy; | |
| 214 otherwise | |
| 215 error('No such boundary: boundary = %s',boundary); | |
| 216 end | |
| 217 end | |
| 218 | |
| 219 function N = size(obj) | |
| 220 N = prod(obj.m); | |
| 221 end | |
| 222 | |
| 223 end | |
| 224 | |
| 225 methods(Static) | |
| 226 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
| 227 % and bound_v of scheme schm_v. | |
| 228 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
| 229 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 230 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 231 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 232 end | |
| 233 end | |
| 234 end |