Mercurial > repos > public > sbplib

comparison +scheme/Burgers1d.m @ 1148:0a5503a08a36 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement the correct Dirichlet conditions for Burgers1d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 24 Jan 2019 09:01:10 +0100
parents 8537fdd6830a
children 635386c073b9
comparison
equal deleted inserted replaced
1147:c322a9454d75 1148:0a5503a08a36
28 ops = opSet(m, xlim, order); 28 ops = opSet(m, xlim, order);
29 29
30 if (isequal(opSet, @sbp.D1Upwind)) 30 if (isequal(opSet, @sbp.D1Upwind))
31 obj.D1 = (ops.Dp + ops.Dm)/2; 31 obj.D1 = (ops.Dp + ops.Dm)/2;
32 DissOp = (ops.Dm - ops.Dp)/2; 32 DissOp = (ops.Dm - ops.Dp)/2;
33 switch pde_form
34 case 'skew-symmetric'
35 obj.D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1 + fluxSplitting(v)*DissOp)*v);
36 case 'conservative'
37 obj.D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
38 end
33 else 39 else
34 obj.D1 = ops.D1; 40 obj.D1 = ops.D1;
41 switch pde_form
42 case 'skew-symmetric'
43 obj.D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
44 case 'conservative'
45 obj.D = @(v) -1/2*obj.D1*v.*v;
46 end
35 end 47 end
36
37 switch pde_form
38 case 'skew-symmetric'
39 if (isequal(opSet, @sbp.D1Upwind))
40 D = @(v) -(1/3*obj.D1*v.*v + (1/3*spdiag(v)*obj.D1 + fluxSplitting(v)*DissOp)*v);
41 else
42 D = @(v) -(1/3*obj.D1*v.*v + 1/3*spdiag(v)*obj.D1*v);
43 end
44 case 'conservative'
45 if (isequal(opSet, @sbp.D1Upwind))
46 D = @(v) -(1/2*obj.D1*v.*v + fluxSplitting(v)*DissOp*v);
47 else
48 D = @(v) -(1/2*obj.D1*v.*v);
49 end
50 otherwise
51 error('Not supported', pde_form);
52 end
53
54 obj.grid = g; 48 obj.grid = g;
55 49
56 obj.D = D;
57 obj.H = ops.H; 50 obj.H = ops.H;
58 obj.Hi = ops.HI; 51 obj.Hi = ops.HI;
59 52
60 obj.e_l = ops.e_l; 53 obj.e_l = ops.e_l;
61 obj.e_r = ops.e_r; 54 obj.e_r = ops.e_r;
68 % Closure functions return the operators applied to the own doamin to close the boundary 61 % Closure functions return the operators applied to the own doamin to close the boundary
69 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. 62 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
70 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 63 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
71 % type is a string specifying the type of boundary condition if there are several. 64 % type is a string specifying the type of boundary condition if there are several.
72 function [closure, penalty] = boundary_condition(obj, boundary, type) 65 function [closure, penalty] = boundary_condition(obj, boundary, type)
73 default_arg('type','robin'); 66 default_arg('type','dirichlet');
74 [e, index, s] = obj.get_boundary_ops(boundary); 67 [e, index, s] = obj.get_boundary_ops(boundary);
75 switch type 68 switch type
76 % Stable dirhclet-like boundary conditions ((u+-abs(u))*u/3)) with +- at left/right boundary 69 % Stable dirchlet-like boundary conditions (u+-abs(u))*u/3
70 % with +- at left/right boundaries
77 case {'D', 'd', 'dirichlet', 'Dirichlet'} 71 case {'D', 'd', 'dirichlet', 'Dirichlet'}
78 tau = s*e; 72 % tau = s*e;
79 closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index)); 73 % closure = @(v) obj.Hi*tau*(((v(index)-s*abs(v(index)))/3)*v(index));
80 penalty = -obj.Hi*tau; 74 % penalty = -obj.Hi*tau;
75
76 magnitude = 2/3;
77 tau = @(v) s*magnitude*obj.Hi*e*(v(index)-s*abs(v(index)))/2;
78 closure = @(v) tau(v)*v(index);
79 penalty = @(v) -tau(v);
81 otherwise 80 otherwise
82 error('No such boundary condition: type = %s',type); 81 error('No such boundary condition: type = %s',type);
83 end 82 end
84 end 83 end
85 84