Mercurial > repos > public > sbplib

comparison +scheme/Laplace1d.m @ 1033:037f203b9bf5 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with branch feature/advectioRV to utilize the +rv package
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:44:12 +0100
parents cab047de7f5d
children 2b1b944deae1 c12b84fe9b00
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
1 classdef Laplace1d < scheme.Scheme
2 properties
3 grid
4 order % Order accuracy for the approximation
5
6 D % non-stabalized scheme operator
7 H % Discrete norm
8 M % Derivative norm
9 a
10
11 D2
12 Hi
13 e_l
14 e_r
15 d_l
16 d_r
17 gamm
18 end
19
20 methods
21 function obj = Laplace1d(grid, order, a)
22 default_arg('a', 1);
23
24 assertType(grid, 'grid.Cartesian');
25
26 ops = sbp.D2Standard(grid.size(), grid.lim{1}, order);
27
28 obj.D2 = sparse(ops.D2);
29 obj.H = sparse(ops.H);
30 obj.Hi = sparse(ops.HI);
31 obj.M = sparse(ops.M);
32 obj.e_l = sparse(ops.e_l);
33 obj.e_r = sparse(ops.e_r);
34 obj.d_l = -sparse(ops.d1_l);
35 obj.d_r = sparse(ops.d1_r);
36
37
38 obj.grid = grid;
39 obj.order = order;
40
41 obj.a = a;
42 obj.D = a*obj.D2;
43
44 obj.gamm = grid.h*ops.borrowing.M.S;
45 end
46
47
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.
50 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
51 % type is a string specifying the type of boundary condition if there are several.
52 % data is a function returning the data that should be applied at the boundary.
53 % neighbour_scheme is an instance of Scheme that should be interfaced to.
54 % neighbour_boundary is a string specifying which boundary to interface to.
55 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
56 default_arg('type','neumann');
57 default_arg('data',0);
58
59 [e,d,s] = obj.get_boundary_ops(boundary);
60
61 switch type
62 % Dirichlet boundary condition
63 case {'D','dirichlet'}
64 tuning = 1.1;
65 tau1 = -tuning/obj.gamm;
66 tau2 = 1;
67
68 tau = tau1*e + tau2*d;
69
70 closure = obj.a*obj.Hi*tau*e';
71 penalty = obj.a*obj.Hi*tau;
72
73 % Neumann boundary condition
74 case {'N','neumann'}
75 tau = -e;
76
77 closure = obj.a*obj.Hi*tau*d';
78 penalty = -obj.a*obj.Hi*tau;
79
80 % Unknown, boundary condition
81 otherwise
82 error('No such boundary condition: type = %s',type);
83 end
84 end
85
86 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
87 % u denotes the solution in the own domain
88 % v denotes the solution in the neighbour domain
89
90 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
91 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
92
93
94 a_u = obj.a;
95 a_v = neighbour_scheme.a;
96
97 gamm_u = obj.gamm;
98 gamm_v = neighbour_scheme.gamm;
99
100 tuning = 1.1;
101
102 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning;
103 tau2 = 1/2*a_u;
104 sig1 = -1/2;
105 sig2 = 0;
106
107 tau = tau1*e_u + tau2*d_u;
108 sig = sig1*e_u + sig2*d_u;
109
110 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u');
111 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v');
112 end
113
114 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
115 % The right boundary is considered the positive boundary
116 function [e,d,s] = get_boundary_ops(obj,boundary)
117 switch boundary
118 case 'l'
119 e = obj.e_l;
120 d = obj.d_l;
121 s = -1;
122 case 'r'
123 e = obj.e_r;
124 d = obj.d_r;
125 s = 1;
126 otherwise
127 error('No such boundary: boundary = %s',boundary);
128 end
129 end
130
131 function N = size(obj)
132 N = obj.grid.size();
133 end
134
135 end
136
137 methods(Static)
138 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
139 % and bound_v of scheme schm_v.
140 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
141 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
142 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
143 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
144 end
145 end
146 end