Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 814:3a5e635a93fd feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme for 1D Burgers equation - Add scheme for discretizing the 1D burgers equation, with spatially variable coefficients.
author Vidar Stiernstrom <vidar.stiernstrom@it.uu.se>
date Mon, 03 Sep 2018 14:50:27 +0200
parents
children fae41958af4f
comparison
equal deleted inserted replaced
580:2ce903f28193 814:3a5e635a93fd
1 classdef Burgers1D < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 x % Grid
6 order % Order accuracy for the approximation
7
8 D % Non-stabalized scheme operator
9 M % Derivative norm
10 H % Discrete norm
11 Hi % Norm inverse
12 e_l
13 e_r
14 d_l
15 d_r
16 params % Parameters for the coefficient matrices
17 end
18
19 methods
20 function obj = Burgers1D(m, order, xlim, params)
21 [x, h] = util.get_grid(xlim{:},m);
22 ops = sbp.D2Variable(m, xlim, order);
23
24 obj.m = m;
25 obj.h = h;
26 obj.order = order;
27 obj.x = x;
28
29 D1 = ops.D1;
30 obj.D = @(v)(-1/3*v.*D1*v - 1/3*D1*v.^2 + ops.D2(params.eps)*v);
31 obj.M = ops.M;
32 obj.H = ops.H;
33 obj.Hi = ops.HI;
34 obj.e_l = ops.e_l;
35 obj.e_r = ops.e_r;
36 obj.d_l = ops.d1_l;
37 obj.d_r = ops.d1_r;
38 obj.params = params;
39 end
40
41 % Closure functions return the opertors applied to the own doamin to close the boundary
42 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other domain.
43 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
44 % type is a string specifying the type of boundary condition if there are several.
45 % data is a function returning the data that should be applied at the boundary.
46 % neighbour_scheme is an instance of Scheme that should be interfaced to.
47 % neighbour_boundary is a string specifying which boundary to interface to.
48 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
49 default_arg('type','robin');
50 default_arg('data',0);
51 [e, s, d] = obj.get_boundary_ops(boundary);
52 switch type
53 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
54 case {'R','robin'}
55 p = s*obj.Hi*e;
56 closure = @(v) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*obj.params.eps*d'*v);
57 switch class(data)
58 case 'double'
59 penalty = s*p*data;
60 case 'function_handle'
61 penalty = @(t) s*p*data(t);
62 otherwise
63 error('Wierd data argument!')
64 end
65 % Unknown, boundary condition
66 otherwise
67 error('No such boundary condition: type = %s',type);
68 end
69 end
70
71 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
72 % The right boundary is considered the positive boundary
73 function [e, s, d] = get_boundary_ops(obj,boundary)
74 switch boundary
75 case 'l'
76 e = obj.e_l;
77 s = -1;
78 d = obj.d_l;
79 case 'r'
80 e = obj.e_r;
81 s = 1;
82 d = obj.d_r;
83 otherwise
84 error('No such boundary: boundary = %s',boundary);
85 end
86 end
87
88 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
89 error('An interface function does not exist yet');
90 end
91
92 function N = size(obj)
93 N = obj.m;
94 end
95
96 end
97 end