Mercurial > repos > public > sbplib
annotate +scheme/Beam2d.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 | a35ed1d124d3 |
| children | 78db023a7fe3 |
| rev | line source |
|---|---|
|
141
cb2b12246b7e
Fixed path to some superclasses.
Jonatan Werpers <jonatan@werpers.com>
parents:
125
diff
changeset
|
1 classdef Beam2d < scheme.Scheme |
| 0 | 2 properties |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
3 grid |
| 0 | 4 order % Order accuracy for the approximation |
| 5 | |
| 6 D % non-stabalized scheme operator | |
| 7 M % Derivative norm | |
| 8 alpha | |
| 9 | |
| 10 H % Discrete norm | |
| 11 Hi | |
| 12 H_x, H_y % Norms in the x and y directions | |
| 13 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 14 Hi_x, Hi_y | |
| 15 Hix, Hiy | |
| 16 e_w, e_e, e_s, e_n | |
| 17 d1_w, d1_e, d1_s, d1_n | |
| 18 d2_w, d2_e, d2_s, d2_n | |
| 19 d3_w, d3_e, d3_s, d3_n | |
| 20 gamm_x, gamm_y | |
| 21 delt_x, delt_y | |
| 22 end | |
| 23 | |
| 24 methods | |
|
125
d52e5cdb6eff
Fixed some class name, file name errors.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
25 function obj = Beam2d(m,lim,order,alpha,opsGen) |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
26 default_arg('alpha',1); |
| 0 | 27 default_arg('opsGen',@sbp.Higher); |
| 28 | |
|
176
d095b5396103
Fixed some bugs in Beam schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
175
diff
changeset
|
29 if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 2 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
30 error('Grid must be 2d cartesian'); |
| 0 | 31 end |
| 32 | |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
33 obj.grid = grid; |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
34 obj.alpha = alpha; |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
35 obj.order = order; |
| 0 | 36 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
37 m_x = grid.m(1); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
38 m_y = grid.m(2); |
| 0 | 39 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
40 h = grid.scaling(); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
41 h_x = h(1); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
42 h_y = h(2); |
| 0 | 43 |
| 44 ops_x = opsGen(m_x,h_x,order); | |
| 45 ops_y = opsGen(m_y,h_y,order); | |
| 46 | |
| 47 I_x = speye(m_x); | |
| 48 I_y = speye(m_y); | |
| 49 | |
| 50 D4_x = sparse(ops_x.derivatives.D4); | |
| 51 H_x = sparse(ops_x.norms.H); | |
| 52 Hi_x = sparse(ops_x.norms.HI); | |
| 53 e_l_x = sparse(ops_x.boundary.e_1); | |
| 54 e_r_x = sparse(ops_x.boundary.e_m); | |
| 55 d1_l_x = sparse(ops_x.boundary.S_1); | |
| 56 d1_r_x = sparse(ops_x.boundary.S_m); | |
| 57 d2_l_x = sparse(ops_x.boundary.S2_1); | |
| 58 d2_r_x = sparse(ops_x.boundary.S2_m); | |
| 59 d3_l_x = sparse(ops_x.boundary.S3_1); | |
| 60 d3_r_x = sparse(ops_x.boundary.S3_m); | |
| 61 | |
| 62 D4_y = sparse(ops_y.derivatives.D4); | |
| 63 H_y = sparse(ops_y.norms.H); | |
| 64 Hi_y = sparse(ops_y.norms.HI); | |
| 65 e_l_y = sparse(ops_y.boundary.e_1); | |
| 66 e_r_y = sparse(ops_y.boundary.e_m); | |
| 67 d1_l_y = sparse(ops_y.boundary.S_1); | |
| 68 d1_r_y = sparse(ops_y.boundary.S_m); | |
| 69 d2_l_y = sparse(ops_y.boundary.S2_1); | |
| 70 d2_r_y = sparse(ops_y.boundary.S2_m); | |
| 71 d3_l_y = sparse(ops_y.boundary.S3_1); | |
| 72 d3_r_y = sparse(ops_y.boundary.S3_m); | |
| 73 | |
| 74 | |
| 75 D4 = kr(D4_x, I_y) + kr(I_x, D4_y); | |
| 76 | |
| 77 % Norms | |
| 78 obj.H = kr(H_x,H_y); | |
| 79 obj.Hx = kr(H_x,I_x); | |
| 80 obj.Hy = kr(I_x,H_y); | |
| 81 obj.Hix = kr(Hi_x,I_y); | |
| 82 obj.Hiy = kr(I_x,Hi_y); | |
| 83 obj.Hi = kr(Hi_x,Hi_y); | |
| 84 | |
| 85 % Boundary operators | |
| 86 obj.e_w = kr(e_l_x,I_y); | |
| 87 obj.e_e = kr(e_r_x,I_y); | |
| 88 obj.e_s = kr(I_x,e_l_y); | |
| 89 obj.e_n = kr(I_x,e_r_y); | |
| 90 obj.d1_w = kr(d1_l_x,I_y); | |
| 91 obj.d1_e = kr(d1_r_x,I_y); | |
| 92 obj.d1_s = kr(I_x,d1_l_y); | |
| 93 obj.d1_n = kr(I_x,d1_r_y); | |
| 94 obj.d2_w = kr(d2_l_x,I_y); | |
| 95 obj.d2_e = kr(d2_r_x,I_y); | |
| 96 obj.d2_s = kr(I_x,d2_l_y); | |
| 97 obj.d2_n = kr(I_x,d2_r_y); | |
| 98 obj.d3_w = kr(d3_l_x,I_y); | |
| 99 obj.d3_e = kr(d3_r_x,I_y); | |
| 100 obj.d3_s = kr(I_x,d3_l_y); | |
| 101 obj.d3_n = kr(I_x,d3_r_y); | |
| 102 | |
| 103 obj.D = alpha*D4; | |
| 104 | |
| 105 obj.gamm_x = h_x*ops_x.borrowing.N.S2/2; | |
| 106 obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2; | |
| 107 | |
| 108 obj.gamm_y = h_y*ops_y.borrowing.N.S2/2; | |
| 109 obj.delt_y = h_y^3*ops_y.borrowing.N.S3/2; | |
| 110 end | |
| 111 | |
| 112 | |
| 113 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 114 % 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 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 116 % type is a string specifying the type of boundary condition if there are several. | |
| 117 % data is a function returning the data that should be applied at the boundary. | |
| 118 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 119 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 120 function [closure, penalty_e,penalty_d] = boundary_condition(obj,boundary,type,data) | |
| 121 default_arg('type','dn'); | |
| 122 default_arg('data',0); | |
| 123 | |
| 124 [e,d1,d2,d3,s,gamm,delt,halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 125 | |
| 126 switch type | |
| 127 % Dirichlet-neumann boundary condition | |
| 128 case {'dn'} | |
| 129 alpha = obj.alpha; | |
| 130 | |
| 131 % tau1 < -alpha^2/gamma | |
| 132 tuning = 1.1; | |
| 133 | |
| 134 tau1 = tuning * alpha/delt; | |
| 135 tau4 = s*alpha; | |
| 136 | |
| 137 sig2 = tuning * alpha/gamm; | |
| 138 sig3 = -s*alpha; | |
| 139 | |
| 140 tau = tau1*e+tau4*d3; | |
| 141 sig = sig2*d1+sig3*d2; | |
| 142 | |
| 143 closure = halfnorm_inv*(tau*e' + sig*d1'); | |
| 144 | |
| 145 pp_e = halfnorm_inv*tau; | |
| 146 pp_d = halfnorm_inv*sig; | |
| 147 switch class(data) | |
| 148 case 'double' | |
| 149 penalty_e = pp_e*data; | |
| 150 penalty_d = pp_d*data; | |
| 151 case 'function_handle' | |
| 152 penalty_e = @(t)pp_e*data(t); | |
| 153 penalty_d = @(t)pp_d*data(t); | |
| 154 otherwise | |
| 155 error('Wierd data argument!') | |
| 156 end | |
| 157 | |
| 158 % Unknown, boundary condition | |
| 159 otherwise | |
| 160 error('No such boundary condition: type = %s',type); | |
| 161 end | |
| 162 end | |
| 163 | |
|
944
a35ed1d124d3
Change from opts to type in a few schemes
Jonatan Werpers <jonatan@werpers.com>
parents:
910
diff
changeset
|
164 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) |
| 0 | 165 % u denotes the solution in the own domain |
| 166 % v denotes the solution in the neighbour domain | |
| 167 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 168 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 169 | |
| 170 tuning = 2; | |
| 171 | |
| 172 alpha_u = obj.alpha; | |
| 173 alpha_v = neighbour_scheme.alpha; | |
| 174 | |
| 175 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; | |
| 176 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; | |
| 177 tau4 = s_u*alpha_u/2; | |
| 178 | |
| 179 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning; | |
| 180 sig3 = -s_u*alpha_u/2; | |
| 181 | |
| 182 phi2 = s_u*1/2; | |
| 183 | |
| 184 psi1 = -s_u*1/2; | |
| 185 | |
| 186 tau = tau1*e_u + tau4*d3_u; | |
| 187 sig = sig2*d1_u + sig3*d2_u ; | |
| 188 phi = phi2*d1_u ; | |
| 189 psi = psi1*e_u ; | |
| 190 | |
| 191 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); | |
| 192 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); | |
| 193 end | |
| 194 | |
| 195 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 196 % The right boundary is considered the positive boundary | |
| 197 function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) | |
| 198 switch boundary | |
| 199 case 'w' | |
| 200 e = obj.e_w; | |
| 201 d1 = obj.d1_w; | |
| 202 d2 = obj.d2_w; | |
| 203 d3 = obj.d3_w; | |
| 204 s = -1; | |
| 205 gamm = obj.gamm_x; | |
| 206 delt = obj.delt_x; | |
| 207 halfnorm_inv = obj.Hix; | |
| 208 case 'e' | |
| 209 e = obj.e_e; | |
| 210 d1 = obj.d1_e; | |
| 211 d2 = obj.d2_e; | |
| 212 d3 = obj.d3_e; | |
| 213 s = 1; | |
| 214 gamm = obj.gamm_x; | |
| 215 delt = obj.delt_x; | |
| 216 halfnorm_inv = obj.Hix; | |
| 217 case 's' | |
| 218 e = obj.e_s; | |
| 219 d1 = obj.d1_s; | |
| 220 d2 = obj.d2_s; | |
| 221 d3 = obj.d3_s; | |
| 222 s = -1; | |
| 223 gamm = obj.gamm_y; | |
| 224 delt = obj.delt_y; | |
| 225 halfnorm_inv = obj.Hiy; | |
| 226 case 'n' | |
| 227 e = obj.e_n; | |
| 228 d1 = obj.d1_n; | |
| 229 d2 = obj.d2_n; | |
| 230 d3 = obj.d3_n; | |
| 231 s = 1; | |
| 232 gamm = obj.gamm_y; | |
| 233 delt = obj.delt_y; | |
| 234 halfnorm_inv = obj.Hiy; | |
| 235 otherwise | |
| 236 error('No such boundary: boundary = %s',boundary); | |
| 237 end | |
| 238 end | |
| 239 | |
| 240 function N = size(obj) | |
| 241 N = prod(obj.m); | |
| 242 end | |
| 243 | |
| 244 end | |
| 245 end |