Mercurial > repos > public > sbplib

comparison +scheme/Laplace1D.m @ 898:bd79326ebcd0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Mordernize Laplace1d
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 22 Nov 2018 10:44:11 +0100
parents 09c5fbc783d3
children 21394c78c72e
comparison
equal deleted inserted replaced
897:ba7e442ea639 898:bd79326ebcd0
1 classdef Laplace1D < scheme.Scheme 1 classdef Laplace1D < scheme.Scheme
2 properties 2 properties
3 g 3 grid
4 order % Order accuracy for the approximation 4 order % Order accuracy for the approximation
5 5
6 D % non-stabalized scheme operator 6 D % non-stabalized scheme operator
7 H % Discrete norm 7 H % Discrete norm
8 M % Derivative norm 8 M % Derivative norm
16 d_r 16 d_r
17 gamm 17 gamm
18 end 18 end
19 19
20 methods 20 methods
21 function obj = Laplace1D(g, order, a) 21 function obj = Laplace1D(grid, order, a)
22 default_arg('a', 1); 22 default_arg('a', 1);
23 23
24 assertType(g, 'grid.Cartesian'); 24 assertType(grid, 'grid.Cartesian');
25 25
26 ops = sbp.Ordinary(g.size(), g.h, order); 26 ops = sbp.D2Standard(grid.size(), grid.lim{1}, order);
27 27
28 obj.D2 = sparse(ops.derivatives.D2); 28 obj.D2 = sparse(ops.D2);
29 obj.H = sparse(ops.norms.H); 29 obj.H = sparse(ops.H);
30 obj.Hi = sparse(ops.norms.HI); 30 obj.Hi = sparse(ops.HI);
31 obj.M = sparse(ops.norms.M); 31 obj.M = sparse(ops.M);
32 obj.e_l = sparse(ops.boundary.e_1); 32 obj.e_l = sparse(ops.e_l);
33 obj.e_r = sparse(ops.boundary.e_m); 33 obj.e_r = sparse(ops.e_r);
34 obj.d_l = sparse(ops.boundary.S_1); 34 obj.d_l = -sparse(ops.d1_l);
35 obj.d_r = sparse(ops.boundary.S_m); 35 obj.d_r = sparse(ops.d1_r);
36 36
37 37
38 obj.g = g; 38 obj.grid = grid;
39 obj.order = order; 39 obj.order = order;
40 40
41 obj.a = a; 41 obj.a = a;
42 obj.D = a*obj.D2; 42 obj.D = a*obj.D2;
43 43
44 obj.gamm = h*ops.borrowing.M.S; 44 obj.gamm = grid.h*ops.borrowing.M.S;
45 end 45 end
46 46
47 47
48 % Closure functions return the opertors applied to the own doamin to close the boundary 48 % Closure functions return the opertors applied to the own doamin to close the boundary
49 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 49 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
59 [e,d,s] = obj.get_boundary_ops(boundary); 59 [e,d,s] = obj.get_boundary_ops(boundary);
60 60
61 switch type 61 switch type
62 % Dirichlet boundary condition 62 % Dirichlet boundary condition
63 case {'D','dirichlet'} 63 case {'D','dirichlet'}
64 alpha = obj.alpha; 64 tuning = 1.1;
65 tau1 = -tuning/obj.gamm;
66 tau2 = 1;
65 67
66 % tau1 < -alpha^2/gamma 68 tau = tau1*e + tau2*d;
67 tuning = 1.1;
68 tau1 = -tuning*alpha/obj.gamm;
69 tau2 = s*alpha;
70 69
71 p = tau1*e + tau2*d; 70 closure = obj.a*obj.Hi*tau*e';
72 71 penalty = obj.a*obj.Hi*tau;
73 closure = obj.Hi*p*e';
74
75 pp = obj.Hi*p;
76 switch class(data)
77 case 'double'
78 penalty = pp*data;
79 case 'function_handle'
80 penalty = @(t)pp*data(t);
81 otherwise
82 error('Wierd data argument!')
83 end
84
85 72
86 % Neumann boundary condition 73 % Neumann boundary condition
87 case {'N','neumann'} 74 case {'N','neumann'}
88 alpha = obj.alpha; 75 tau = -e;
89 tau1 = -s*alpha;
90 tau2 = 0;
91 tau = tau1*e + tau2*d;
92 76
93 closure = obj.Hi*tau*d'; 77 closure = obj.a*obj.Hi*tau*d';
94 78 penalty = -obj.a*obj.Hi*tau;
95 pp = obj.Hi*tau;
96 switch class(data)
97 case 'double'
98 penalty = pp*data;
99 case 'function_handle'
100 penalty = @(t)pp*data(t);
101 otherwise
102 error('Wierd data argument!')
103 end
104 79
105 % Unknown, boundary condition 80 % Unknown, boundary condition
106 otherwise 81 otherwise
107 error('No such boundary condition: type = %s',type); 82 error('No such boundary condition: type = %s',type);
108 end 83 end
109 end 84 end
110 85
111 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 86 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
112 % u denotes the solution in the own domain 87 % u denotes the solution in the own domain
113 % v denotes the solution in the neighbour domain 88 % v denotes the solution in the neighbour domain
89
114 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); 90 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
115 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 91 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
116 92
117 tuning = 1.1;
118 93
119 alpha_u = obj.alpha; 94 a_u = obj.a;
120 alpha_v = neighbour_scheme.alpha; 95 a_v = neighbour_scheme.a;
121 96
122 gamm_u = obj.gamm; 97 gamm_u = obj.gamm;
123 gamm_v = neighbour_scheme.gamm; 98 gamm_v = neighbour_scheme.gamm;
124 99
125 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v) 100 tuning = 1.1;
126 101
127 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning; 102 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning;
128 tau2 = s_u*1/2*alpha_u; 103 tau2 = 1/2*a_u;
129 sig1 = s_u*(-1/2); 104 sig1 = -1/2;
130 sig2 = 0; 105 sig2 = 0;
131 106
132 tau = tau1*e_u + tau2*d_u; 107 tau = tau1*e_u + tau2*d_u;
133 sig = sig1*e_u + sig2*d_u; 108 sig = sig1*e_u + sig2*d_u;
134 109
135 closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u'); 110 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
136 penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v'); 111 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
137 end 112 end
138 113
139 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 114 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
140 % The right boundary is considered the positive boundary 115 % The right boundary is considered the positive boundary
141 function [e,d,s] = get_boundary_ops(obj,boundary) 116 function [e,d,s] = get_boundary_ops(obj,boundary)
152 error('No such boundary: boundary = %s',boundary); 127 error('No such boundary: boundary = %s',boundary);
153 end 128 end
154 end 129 end
155 130
156 function N = size(obj) 131 function N = size(obj)
157 N = obj.m; 132 N = obj.grid.size();
158 end 133 end
159 134
160 end 135 end
161 136
162 methods(Static) 137 methods(Static)