Mercurial > repos > public > sbplib

comparison +scheme/Utux.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 d6ab5ceba496
children 8a9393084b30
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
1 classdef Utux < scheme.Scheme 1 classdef Utux < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 x % Grid 5 grid % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7
8 a % Wave speed
9 % Can either be a constant or function handle.
7 10
8 H % Discrete norm 11 H % Discrete norm
9 D 12 D
10 13
11 D1 14 D1
12 Hi 15 Hi
13 e_l 16 e_l
14 e_r 17 e_r
15 v0
16 end 18 end
17 19
18 20
19 methods 21 methods
20 function obj = Utux(m,xlim,order,operator) 22 function obj = Utux(g, order, opSet, a, fluxSplitting)
21 default_arg('a',1); 23 default_arg('opSet',@sbp.D2Standard);
22 24 default_arg('a',1);
23 %Old operators 25 default_arg('fluxSplitting',[]);
24 % [x, h] = util.get_grid(xlim{:},m);
25 %ops = sbp.Ordinary(m,h,order);
26
27
28 switch operator
29 case 'NonEquidistant'
30 ops = sbp.D1Nonequidistant(m,xlim,order);
31 obj.D1 = ops.D1;
32 case 'Standard'
33 ops = sbp.D2Standard(m,xlim,order);
34 obj.D1 = ops.D1;
35 case 'Upwind'
36 ops = sbp.D1Upwind(m,xlim,order);
37 obj.D1 = ops.Dm;
38 otherwise
39 error('Unvalid operator')
40 end
41 obj.x=ops.x;
42 26
43 27 assertType(g, 'grid.Cartesian');
28 if isa(a, 'function_handle')
29 obj.a = spdiag(grid.evalOn(g, a));
30 else
31 obj.a = a;
32 end
33
34 m = g.size();
35 xl = g.getBoundary('l');
36 xr = g.getBoundary('r');
37 xlim = {xl, xr};
38
39 ops = opSet(m, xlim, order);
40
41 if (isequal(opSet, @sbp.D1Upwind))
42 obj.D1 = (ops.Dp + ops.Dm)/2;
43 DissOp = (ops.Dm - ops.Dp)/2;
44 obj.D = -(obj.a*obj.D1 + fluxSplitting*DissOp);
45 else
46 obj.D1 = ops.D1;
47 obj.D = -obj.a*obj.D1;
48 end
49
50 obj.grid = g;
51
44 obj.H = ops.H; 52 obj.H = ops.H;
45 obj.Hi = ops.HI; 53 obj.Hi = ops.HI;
46 54
47 obj.e_l = ops.e_l; 55 obj.e_l = ops.e_l;
48 obj.e_r = ops.e_r; 56 obj.e_r = ops.e_r;
49 obj.D=obj.D1;
50 57
51 obj.m = m; 58 obj.m = m;
52 obj.h = ops.h; 59 obj.h = ops.h;
53 obj.order = order; 60 obj.order = order;
54 obj.x = ops.x;
55
56 end 61 end
57 % Closure functions return the opertors applied to the own doamin to close the boundary 62 % Closure functions return the opertors applied to the own doamin to close the boundary
58 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 63 % 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 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 64 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
60 % type is a string specifying the type of boundary condition if there are several. 65 % type is a string specifying the type of boundary condition if there are several.
61 % data is a function returning the data that should be applied at the boundary. 66 % data is a function returning the data that should be applied at the boundary.
62 % neighbour_scheme is an instance of Scheme that should be interfaced to. 67 % neighbour_scheme is an instance of Scheme that should be interfaced to.
63 % neighbour_boundary is a string specifying which boundary to interface to. 68 % neighbour_boundary is a string specifying which boundary to interface to.
64 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 69 function [closure, penalty] = boundary_condition(obj,boundary,type)
65 default_arg('type','neumann'); 70 default_arg('type','dirichlet');
66 default_arg('data',0); 71 sigma_left = -1; % Scalar penalty parameter for left boundary
67 tau =-1*obj.e_l; 72 sigma_right = 1; % Scalar penalty parameter for right boundary
68 closure = obj.Hi*tau*obj.e_l'; 73 switch boundary
69 penalty = 0*obj.e_l; 74 % Can only specify boundary condition where there is inflow
70 75 % Extract the postivie resp. negative part of a, for the left
76 % resp. right boundary, and set other values of a to zero.
77 % Then the closure will effectively only contribute to inflow boundaries
78 case {'l','L','left','Left'}
79 a_inflow = obj.a;
80 a_inflow(a_inflow < 0) = 0;
81 tau = sigma_left*a_inflow*obj.e_l;
82 closure = obj.Hi*tau*obj.e_l';
83 case {'r','R','right','Right'}
84 a_inflow = obj.a;
85 a_inflow(a_inflow > 0) = 0;
86 tau = sigma_right*a_inflow*obj.e_r;
87 closure = obj.Hi*tau*obj.e_r';
88 end
89 penalty = -obj.Hi*tau;
90
71 end 91 end
72 92
73 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 93 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
74 error('An interface function does not exist yet'); 94 switch boundary
95 % Upwind coupling
96 case {'l','left'}
97 tau = -1*obj.a*obj.e_l;
98 closure = obj.Hi*tau*obj.e_l';
99 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
100 case {'r','right'}
101 tau = 0*obj.a*obj.e_r;
102 closure = obj.Hi*tau*obj.e_r';
103 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
104 end
105
75 end 106 end
76 107
77 function N = size(obj) 108 function N = size(obj)
78 N = obj.m; 109 N = obj.m;
79 end 110 end
80 111
81 end 112 end
82 113
83 methods(Static) 114 methods(Static)
84 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u 115 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
85 % and bound_v of scheme schm_v. 116 % and bound_v of scheme schm_v.
86 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') 117 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
87 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) 118 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
88 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); 119 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
89 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); 120 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
90 end 121 end
91 end 122 end