Mercurial > repos > public > sbplib
annotate +scheme/Laplace1d.m @ 1037:2d7ba44340d0 feature/burgers1d
Pass scheme specific parameters as cell array. This will enabale constructDiffOps to be more general. In addition, allow for schemes returning function handles as diffOps, which is currently how non-linear schemes such as Burgers1d are implemented.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Fri, 18 Jan 2019 09:02:02 +0100 |
| parents | cab047de7f5d |
| children | 2b1b944deae1 c12b84fe9b00 |
| rev | line source |
|---|---|
|
950
cab047de7f5d
Rename *2D schemes to *2d
Jonatan Werpers <jonatan@werpers.com>
parents:
946
diff
changeset
|
1 classdef Laplace1d < scheme.Scheme |
| 0 | 2 properties |
| 898 | 3 grid |
| 0 | 4 order % Order accuracy for the approximation |
| 5 | |
| 6 D % non-stabalized scheme operator | |
| 7 H % Discrete norm | |
| 8 M % Derivative norm | |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
9 a |
| 0 | 10 |
| 11 D2 | |
| 12 Hi | |
| 13 e_l | |
| 14 e_r | |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
15 d_l |
|
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
16 d_r |
| 0 | 17 gamm |
| 18 end | |
| 19 | |
| 20 methods | |
|
950
cab047de7f5d
Rename *2D schemes to *2d
Jonatan Werpers <jonatan@werpers.com>
parents:
946
diff
changeset
|
21 function obj = Laplace1d(grid, order, a) |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
22 default_arg('a', 1); |
| 0 | 23 |
| 898 | 24 assertType(grid, 'grid.Cartesian'); |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
25 |
| 898 | 26 ops = sbp.D2Standard(grid.size(), grid.lim{1}, order); |
| 0 | 27 |
| 898 | 28 obj.D2 = sparse(ops.D2); |
| 29 obj.H = sparse(ops.H); | |
| 30 obj.Hi = sparse(ops.HI); | |
| 31 obj.M = sparse(ops.M); | |
| 32 obj.e_l = sparse(ops.e_l); | |
| 33 obj.e_r = sparse(ops.e_r); | |
| 34 obj.d_l = -sparse(ops.d1_l); | |
| 35 obj.d_r = sparse(ops.d1_r); | |
| 0 | 36 |
| 37 | |
| 898 | 38 obj.grid = grid; |
| 0 | 39 obj.order = order; |
| 40 | |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
41 obj.a = a; |
|
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
42 obj.D = a*obj.D2; |
| 0 | 43 |
| 898 | 44 obj.gamm = grid.h*ops.borrowing.M.S; |
| 0 | 45 end |
| 46 | |
| 47 | |
| 48 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 49 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 50 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 51 % type is a string specifying the type of boundary condition if there are several. | |
| 52 % data is a function returning the data that should be applied at the boundary. | |
| 53 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 54 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 55 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
| 56 default_arg('type','neumann'); | |
| 57 default_arg('data',0); | |
| 58 | |
| 59 [e,d,s] = obj.get_boundary_ops(boundary); | |
| 60 | |
| 61 switch type | |
| 62 % Dirichlet boundary condition | |
| 63 case {'D','dirichlet'} | |
| 64 tuning = 1.1; | |
| 898 | 65 tau1 = -tuning/obj.gamm; |
| 66 tau2 = 1; | |
| 0 | 67 |
| 898 | 68 tau = tau1*e + tau2*d; |
| 0 | 69 |
| 898 | 70 closure = obj.a*obj.Hi*tau*e'; |
| 71 penalty = obj.a*obj.Hi*tau; | |
| 0 | 72 |
| 73 % Neumann boundary condition | |
| 74 case {'N','neumann'} | |
| 898 | 75 tau = -e; |
| 0 | 76 |
| 898 | 77 closure = obj.a*obj.Hi*tau*d'; |
| 78 penalty = -obj.a*obj.Hi*tau; | |
| 0 | 79 |
| 80 % Unknown, boundary condition | |
| 81 otherwise | |
| 82 error('No such boundary condition: type = %s',type); | |
| 83 end | |
| 84 end | |
| 85 | |
|
946
706d1c2b4199
Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents:
943
diff
changeset
|
86 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) |
| 0 | 87 % u denotes the solution in the own domain |
| 88 % v denotes the solution in the neighbour domain | |
| 898 | 89 |
| 0 | 90 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); |
| 91 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 92 | |
| 93 | |
| 898 | 94 a_u = obj.a; |
| 95 a_v = neighbour_scheme.a; | |
| 0 | 96 |
| 97 gamm_u = obj.gamm; | |
| 98 gamm_v = neighbour_scheme.gamm; | |
| 99 | |
| 898 | 100 tuning = 1.1; |
|
946
706d1c2b4199
Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents:
943
diff
changeset
|
101 |
| 898 | 102 tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning; |
| 103 tau2 = 1/2*a_u; | |
| 104 sig1 = -1/2; | |
| 0 | 105 sig2 = 0; |
| 106 | |
| 107 tau = tau1*e_u + tau2*d_u; | |
| 108 sig = sig1*e_u + sig2*d_u; | |
| 109 | |
| 898 | 110 closure = obj.Hi*( tau*e_u' + sig*a_u*d_u'); |
| 111 penalty = obj.Hi*(-tau*e_v' + sig*a_v*d_v'); | |
| 0 | 112 end |
| 113 | |
| 114 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 115 % The right boundary is considered the positive boundary | |
| 116 function [e,d,s] = get_boundary_ops(obj,boundary) | |
| 117 switch boundary | |
| 118 case 'l' | |
| 119 e = obj.e_l; | |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
120 d = obj.d_l; |
| 0 | 121 s = -1; |
| 122 case 'r' | |
| 123 e = obj.e_r; | |
|
896
09c5fbc783d3
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
124 d = obj.d_r; |
| 0 | 125 s = 1; |
| 126 otherwise | |
| 127 error('No such boundary: boundary = %s',boundary); | |
| 128 end | |
| 129 end | |
| 130 | |
| 131 function N = size(obj) | |
| 898 | 132 N = obj.grid.size(); |
| 0 | 133 end |
| 134 | |
| 135 end | |
| 136 | |
| 137 methods(Static) | |
| 138 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
| 139 % and bound_v of scheme schm_v. | |
| 140 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
| 141 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 142 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 143 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 144 end | |
| 145 end | |
| 146 end |