Mercurial > repos > public > sbplib
annotate +scheme/Wave.m @ 774:66eb4a2bbb72 feature/grids
Remove default scaling of the system.
The scaling doens't seem to help actual solutions. One example that fails in the flexural code.
With large timesteps the solutions seems to blow up. One particular example is profilePresentation
on the tdb_presentation_figures branch with k = 0.0005
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 18 Jul 2018 15:42:52 -0700 |
| parents | cb2b12246b7e |
| children | 459eeb99130f |
| rev | line source |
|---|---|
|
141
cb2b12246b7e
Fixed path to some superclasses.
Jonatan Werpers <jonatan@werpers.com>
parents:
125
diff
changeset
|
1 classdef Wave < scheme.Scheme |
| 0 | 2 properties |
| 3 m % Number of points in each direction, possibly a vector | |
| 4 h % Grid spacing | |
| 5 x % Grid | |
| 6 order % Order accuracy for the approximation | |
| 7 | |
| 8 D % non-stabalized scheme operator | |
| 9 H % Discrete norm | |
| 10 M % Derivative norm | |
| 11 alpha | |
| 12 | |
| 13 D2 | |
| 14 Hi | |
| 15 e_l | |
| 16 e_r | |
| 17 d1_l | |
| 18 d1_r | |
| 19 gamm | |
| 20 end | |
| 21 | |
| 22 methods | |
|
125
d52e5cdb6eff
Fixed some class name, file name errors.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
23 function obj = Wave(m,xlim,order,alpha) |
| 0 | 24 default_arg('a',1); |
| 25 [x, h] = util.get_grid(xlim{:},m); | |
| 26 | |
| 27 ops = sbp.Ordinary(m,h,order); | |
| 28 | |
| 29 obj.D2 = sparse(ops.derivatives.D2); | |
| 30 obj.H = sparse(ops.norms.H); | |
| 31 obj.Hi = sparse(ops.norms.HI); | |
| 32 obj.M = sparse(ops.norms.M); | |
| 33 obj.e_l = sparse(ops.boundary.e_1); | |
| 34 obj.e_r = sparse(ops.boundary.e_m); | |
| 35 obj.d1_l = sparse(ops.boundary.S_1); | |
| 36 obj.d1_r = sparse(ops.boundary.S_m); | |
| 37 | |
| 38 | |
| 39 obj.m = m; | |
| 40 obj.h = h; | |
| 41 obj.order = order; | |
| 42 | |
| 43 obj.alpha = alpha; | |
| 44 obj.D = alpha*obj.D2; | |
| 45 obj.x = x; | |
| 46 | |
| 47 obj.gamm = h*ops.borrowing.M.S; | |
| 48 | |
| 49 end | |
| 50 | |
| 51 | |
| 52 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 53 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 54 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 55 % type is a string specifying the type of boundary condition if there are several. | |
| 56 % data is a function returning the data that should be applied at the boundary. | |
| 57 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 58 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 59 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
| 60 default_arg('type','neumann'); | |
| 61 default_arg('data',0); | |
| 62 | |
| 63 [e,d,s] = obj.get_boundary_ops(boundary); | |
| 64 | |
| 65 switch type | |
| 66 % Dirichlet boundary condition | |
| 67 case {'D','dirichlet'} | |
| 68 alpha = obj.alpha; | |
| 69 | |
| 70 % tau1 < -alpha^2/gamma | |
| 71 tuning = 1.1; | |
| 72 tau1 = -tuning*alpha/obj.gamm; | |
| 73 tau2 = s*alpha; | |
| 74 | |
| 75 p = tau1*e + tau2*d; | |
| 76 | |
| 77 closure = obj.Hi*p*e'; | |
| 78 | |
| 79 pp = obj.Hi*p; | |
| 80 switch class(data) | |
| 81 case 'double' | |
| 82 penalty = pp*data; | |
| 83 case 'function_handle' | |
| 84 penalty = @(t)pp*data(t); | |
| 85 otherwise | |
| 86 error('Wierd data argument!') | |
| 87 end | |
| 88 | |
| 89 | |
| 90 % Neumann boundary condition | |
| 91 case {'N','neumann'} | |
| 92 alpha = obj.alpha; | |
| 93 tau1 = -s*alpha; | |
| 94 tau2 = 0; | |
| 95 tau = tau1*e + tau2*d; | |
| 96 | |
| 97 closure = obj.Hi*tau*d'; | |
| 98 | |
| 99 pp = obj.Hi*tau; | |
| 100 switch class(data) | |
| 101 case 'double' | |
| 102 penalty = pp*data; | |
| 103 case 'function_handle' | |
| 104 penalty = @(t)pp*data(t); | |
| 105 otherwise | |
| 106 error('Wierd data argument!') | |
| 107 end | |
| 108 | |
| 109 % Unknown, boundary condition | |
| 110 otherwise | |
| 111 error('No such boundary condition: type = %s',type); | |
| 112 end | |
| 113 end | |
| 114 | |
| 115 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 116 % u denotes the solution in the own domain | |
| 117 % v denotes the solution in the neighbour domain | |
| 118 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); | |
| 119 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 120 | |
| 121 tuning = 1.1; | |
| 122 | |
| 123 alpha_u = obj.alpha; | |
| 124 alpha_v = neighbour_scheme.alpha; | |
| 125 | |
| 126 gamm_u = obj.gamm; | |
| 127 gamm_v = neighbour_scheme.gamm; | |
| 128 | |
| 129 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) | |
| 130 | |
| 131 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; | |
| 132 tau2 = s_u*1/2*alpha_u; | |
| 133 sig1 = s_u*(-1/2); | |
| 134 sig2 = 0; | |
| 135 | |
| 136 tau = tau1*e_u + tau2*d_u; | |
| 137 sig = sig1*e_u + sig2*d_u; | |
| 138 | |
| 139 closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u'); | |
| 140 penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v'); | |
| 141 end | |
| 142 | |
| 143 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 144 % The right boundary is considered the positive boundary | |
| 145 function [e,d,s] = get_boundary_ops(obj,boundary) | |
| 146 switch boundary | |
| 147 case 'l' | |
| 148 e = obj.e_l; | |
| 149 d = obj.d1_l; | |
| 150 s = -1; | |
| 151 case 'r' | |
| 152 e = obj.e_r; | |
| 153 d = obj.d1_r; | |
| 154 s = 1; | |
| 155 otherwise | |
| 156 error('No such boundary: boundary = %s',boundary); | |
| 157 end | |
| 158 end | |
| 159 | |
| 160 function N = size(obj) | |
| 161 N = obj.m; | |
| 162 end | |
| 163 | |
| 164 end | |
| 165 | |
| 166 methods(Static) | |
| 167 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
| 168 % and bound_v of scheme schm_v. | |
| 169 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
| 170 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 171 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 172 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 173 end | |
| 174 end | |
| 175 end |