Mercurial > repos > public > sbplib

comparison +scheme/Burgers1D.m @ 815:fae41958af4f feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for artificial viscosity to the 1d burgers scheme. - Add support for artificial viscosity by parametrizing the discretization operator and boundary closure on the viscosity. - Add a time stepper which evaluates and updates the residual viscosity of the solution.
author Vidar Stiernstrom <vidar.stiernstrom@it.uu.se>
date Thu, 06 Sep 2018 12:43:51 +0200
parents 3a5e635a93fd
children d0934d1143b7
comparison
equal deleted inserted replaced
814:3a5e635a93fd 815:fae41958af4f
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 x % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7 params
7 8
8 D % Non-stabalized scheme operator 9 D % Non-stabalized scheme operator
9 M % Derivative norm 10 M % Derivative norm
10 H % Discrete norm 11 H % Discrete norm
11 Hi % Norm inverse 12 Hi % Norm inverse
12 e_l 13 e_l
13 e_r 14 e_r
14 d_l 15 d_l
15 d_r 16 d_r
16 params % Parameters for the coefficient matrices
17 end 17 end
18 18
19 methods 19 methods
20 function obj = Burgers1D(m, order, xlim, params) 20 function obj = Burgers1D(pde_form, operator_type, order, m, lim, params)
21 [x, h] = util.get_grid(xlim{:},m); 21 [x, h] = util.get_grid(lim{:},m);
22 ops = sbp.D2Variable(m, xlim, order); 22 default_arg('pde_form','skew-symmetric');
23 default_arg('operator_type','narrow');
24
25 switch operator_type
26 case 'narrow'
27 ops = sbp.D2Variable(m, lim, order);
28 D1 = ops.D1;
29 D2 = ops.D2;
30 otherwise
31 error('Other operator types not yet supported', operator_type);
32 end
33
34 switch pde_form
35 case 'skew-symmetric'
36 D = @(v, viscosity) -1/3*v.*D1*v - 1/3*D1*v.^2 + D2(obj.params.eps + viscosity)*v;
37 case 'conservative'
38 D = @(v, viscosity) -1/2*D1*v.^2 + D2(obj.params.eps + viscosity)*v;
39 end
23 40
24 obj.m = m; 41 obj.m = m;
25 obj.h = h; 42 obj.h = h;
26 obj.order = order; 43 obj.order = order;
27 obj.x = x; 44 obj.x = x;
45 obj.params = params;
28 46
29 D1 = ops.D1; 47 %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
30 obj.D = @(v)(-1/3*v.*D1*v - 1/3*D1*v.^2 + ops.D2(params.eps)*v); 48 obj.D = D;
31 obj.M = ops.M; 49 obj.M = ops.M;
32 obj.H = ops.H; 50 obj.H = ops.H;
33 obj.Hi = ops.HI; 51 obj.Hi = ops.HI;
34 obj.e_l = ops.e_l; 52 obj.e_l = ops.e_l;
35 obj.e_r = ops.e_r; 53 obj.e_r = ops.e_r;
36 obj.d_l = ops.d1_l; 54 obj.d_l = ops.d1_l;
37 obj.d_r = ops.d1_r; 55 obj.d_r = ops.d1_r;
38 obj.params = params;
39 end 56 end
40 57
41 % Closure functions return the opertors applied to the own doamin to close the boundary 58 % 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. 59 % 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'. 60 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
51 [e, s, d] = obj.get_boundary_ops(boundary); 68 [e, s, d] = obj.get_boundary_ops(boundary);
52 switch type 69 switch type
53 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary 70 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
54 case {'R','robin'} 71 case {'R','robin'}
55 p = s*obj.Hi*e; 72 p = s*obj.Hi*e;
56 closure = @(v) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*obj.params.eps*d'*v); 73 closure = @(v, viscosity) p*(e'*((v-s*abs(v))/3)*(e'*v) - e'*(obj.params.eps + viscosity)*d'*v);
57 switch class(data) 74 switch class(data)
58 case 'double' 75 case 'double'
59 penalty = s*p*data; 76 penalty = s*p*data;
60 case 'function_handle' 77 case 'function_handle'
61 penalty = @(t) s*p*data(t); 78 penalty = @(t) s*p*data(t);