Mercurial > repos > public > sbplib

comparison +scheme/Utux.m @ 979:7a5e770974ed feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Jan 2019 16:26:00 +0100
parents 6d2167719557
children 2b1b944deae1 d6ab5ceba496 c12b84fe9b00
comparison
equal deleted inserted replaced
933:34b3d092a4d0 979:7a5e770974ed
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 7
8 H % Discrete norm 8 H % Discrete norm
9 D 9 D
10 10
14 e_r 14 e_r
15 v0 15 v0
16 end 16 end
17 17
18 18
19 methods 19 methods
20 function obj = Utux(m,xlim,order,operator) 20 function obj = Utux(g, order, opSet)
21 default_arg('a',1); 21 default_arg('opSet',@sbp.D2Standard);
22
23 %Old operators
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 22
43 23 m = g.size();
24 xl = g.getBoundary('l');
25 xr = g.getBoundary('r');
26 xlim = {xl, xr};
27
28 ops = opSet(m, xlim, order);
29 obj.D1 = ops.D1;
30
31 obj.grid = g;
32
44 obj.H = ops.H; 33 obj.H = ops.H;
45 obj.Hi = ops.HI; 34 obj.Hi = ops.HI;
46 35
47 obj.e_l = ops.e_l; 36 obj.e_l = ops.e_l;
48 obj.e_r = ops.e_r; 37 obj.e_r = ops.e_r;
49 obj.D=obj.D1; 38 obj.D = -obj.D1;
50 39
51 obj.m = m; 40 obj.m = m;
52 obj.h = ops.h; 41 obj.h = ops.h;
53 obj.order = order; 42 obj.order = order;
54 obj.x = ops.x;
55 43
56 end 44 end
57 % Closure functions return the opertors applied to the own doamin to close the boundary 45 % 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. 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.
59 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 47 % 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. 48 % 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. 49 % 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. 50 % neighbour_scheme is an instance of Scheme that should be interfaced to.
63 % neighbour_boundary is a string specifying which boundary to interface to. 51 % neighbour_boundary is a string specifying which boundary to interface to.
64 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 52 function [closure, penalty] = boundary_condition(obj,boundary,type)
65 default_arg('type','neumann'); 53 default_arg('type','dirichlet');
66 default_arg('data',0); 54 tau =-1*obj.e_l;
67 tau =-1*obj.e_l; 55 closure = obj.Hi*tau*obj.e_l';
68 closure = obj.Hi*tau*obj.e_l'; 56 penalty = -obj.Hi*tau;
69 penalty = 0*obj.e_l; 57
70
71 end 58 end
72 59
73 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 60 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
74 error('An interface function does not exist yet'); 61 switch boundary
62 % Upwind coupling
63 case {'l','left'}
64 tau = -1*obj.e_l;
65 closure = obj.Hi*tau*obj.e_l';
66 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
67 case {'r','right'}
68 tau = 0*obj.e_r;
69 closure = obj.Hi*tau*obj.e_r';
70 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
71 end
72
75 end 73 end
76 74
77 function N = size(obj) 75 function N = size(obj)
78 N = obj.m; 76 N = obj.m;
79 end 77 end
80 78
81 end 79 end
82 80
83 methods(Static) 81 methods(Static)
84 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u 82 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
85 % and bound_v of scheme schm_v. 83 % and bound_v of scheme schm_v.
86 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') 84 % [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) 85 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); 86 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
89 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); 87 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
90 end 88 end
91 end 89 end