Mercurial > repos > public > sbplib

comparison +scheme/Utux.m @ 1027:d6ab5ceba496 feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Support variable wave speed and upwind operators in Utux
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 07 Jan 2019 16:35:04 +0100
parents 6d2167719557
children 8a9393084b30
comparison
equal deleted inserted replaced
1026:44c3ea38097e 1027:d6ab5ceba496
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 grid % 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(g, order, opSet) 22 function obj = Utux(g, order, opSet, a, fluxSplitting)
21 default_arg('opSet',@sbp.D2Standard); 23 default_arg('opSet',@sbp.D2Standard);
24 default_arg('a',1);
25 default_arg('fluxSplitting',[]);
26
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
22 33
23 m = g.size(); 34 m = g.size();
24 xl = g.getBoundary('l'); 35 xl = g.getBoundary('l');
25 xr = g.getBoundary('r'); 36 xr = g.getBoundary('r');
26 xlim = {xl, xr}; 37 xlim = {xl, xr};
27 38
28 ops = opSet(m, xlim, order); 39 ops = opSet(m, xlim, order);
29 obj.D1 = ops.D1; 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
30 49
31 obj.grid = g; 50 obj.grid = g;
32 51
33 obj.H = ops.H; 52 obj.H = ops.H;
34 obj.Hi = ops.HI; 53 obj.Hi = ops.HI;
35 54
36 obj.e_l = ops.e_l; 55 obj.e_l = ops.e_l;
37 obj.e_r = ops.e_r; 56 obj.e_r = ops.e_r;
38 obj.D = -obj.D1;
39 57
40 obj.m = m; 58 obj.m = m;
41 obj.h = ops.h; 59 obj.h = ops.h;
42 obj.order = order; 60 obj.order = order;
43
44 end 61 end
45 % 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
46 % 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.
47 % 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'.
48 % 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.
49 % 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.
50 % 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.
51 % neighbour_boundary is a string specifying which boundary to interface to. 68 % neighbour_boundary is a string specifying which boundary to interface to.
52 function [closure, penalty] = boundary_condition(obj,boundary,type) 69 function [closure, penalty] = boundary_condition(obj,boundary,type)
53 default_arg('type','dirichlet'); 70 default_arg('type','dirichlet');
54 tau =-1*obj.e_l; 71 sigma_left = -1; % Scalar penalty parameter for left boundary
55 closure = obj.Hi*tau*obj.e_l'; 72 sigma_right = 1; % Scalar penalty parameter for right boundary
73 switch boundary
74 % Can only specify boundary condition where there is inflow
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
56 penalty = -obj.Hi*tau; 89 penalty = -obj.Hi*tau;
57 90
58 end 91 end
59 92
60 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) 93 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
61 switch boundary 94 switch boundary
62 % Upwind coupling 95 % Upwind coupling
63 case {'l','left'} 96 case {'l','left'}
64 tau = -1*obj.e_l; 97 tau = -1*obj.a*obj.e_l;
65 closure = obj.Hi*tau*obj.e_l'; 98 closure = obj.Hi*tau*obj.e_l';
66 penalty = -obj.Hi*tau*neighbour_scheme.e_r'; 99 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
67 case {'r','right'} 100 case {'r','right'}
68 tau = 0*obj.e_r; 101 tau = 0*obj.a*obj.e_r;
69 closure = obj.Hi*tau*obj.e_r'; 102 closure = obj.Hi*tau*obj.e_r';
70 penalty = -obj.Hi*tau*neighbour_scheme.e_l'; 103 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
71 end 104 end
72 105
73 end 106 end