Mercurial > repos > public > sbplib
annotate +scheme/Schrodinger.m @ 1031:2ef20d00b386 feature/advectionRV
For easier comparison, return both the first order and residual viscosity when evaluating the residual. Add the first order and residual viscosity to the state of the RungekuttaRV time steppers
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 17 Jan 2019 10:25:06 +0100 |
| parents | 706d1c2b4199 |
| children | 337c4d1dcef5 c12b84fe9b00 |
| rev | line source |
|---|---|
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 classdef Schrodinger < scheme.Scheme |
| 0 | 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 H % Discrete norm | |
| 10 M % Derivative norm | |
| 11 alpha | |
| 12 | |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
13 D2 |
| 0 | 14 Hi |
| 15 e_l | |
| 16 e_r | |
| 17 d1_l | |
| 18 d1_r | |
| 19 gamm | |
| 20 end | |
| 21 | |
| 22 methods | |
| 23 % Solving SE in the form u_t = i*u_xx -i*V; | |
| 24 function obj = Schrodinger(m,xlim,order,V) | |
| 25 default_arg('V',0); | |
| 26 | |
| 27 [x, h] = util.get_grid(xlim{:},m); | |
| 28 | |
| 29 ops = sbp.Ordinary(m,h,order); | |
| 30 | |
| 31 obj.D2 = sparse(ops.derivatives.D2); | |
| 32 obj.H = sparse(ops.norms.H); | |
| 33 obj.Hi = sparse(ops.norms.HI); | |
| 34 obj.M = sparse(ops.norms.M); | |
| 35 obj.e_l = sparse(ops.boundary.e_1); | |
| 36 obj.e_r = sparse(ops.boundary.e_m); | |
| 37 obj.d1_l = sparse(ops.boundary.S_1); | |
| 38 obj.d1_r = sparse(ops.boundary.S_m); | |
| 39 | |
| 40 | |
| 41 if isa(V,'function_handle') | |
| 42 V_vec = V(x); | |
| 43 else | |
| 44 V_vec = x*0 + V; | |
| 45 end | |
| 46 | |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
47 V_mat = spdiags(V_vec,0,m,m); |
| 0 | 48 |
|
67
446d67a49cd8
Fixed some errors in scheme.Schrodinger.m
Jonatan Werpers <jonatan@werpers.com>
parents:
65
diff
changeset
|
49 obj.D = 1i * obj.D2 - 1i * V_mat; |
| 0 | 50 |
| 51 obj.m = m; | |
| 52 obj.h = h; | |
| 53 obj.order = order; | |
| 54 | |
| 55 obj.x = x; | |
| 56 end | |
| 57 | |
| 58 | |
| 59 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 60 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 61 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 62 % type is a string specifying the type of boundary condition if there are several. | |
| 63 % data is a function returning the data that should be applied at the boundary. | |
| 64 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 65 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 66 function [closure, penalty] = boundary_condition(obj,boundary,type,data) | |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
67 default_arg('type','dirichlet'); |
| 0 | 68 default_arg('data',0); |
| 69 | |
| 70 [e,d,s] = obj.get_boundary_ops(boundary); | |
| 71 | |
| 72 switch type | |
| 73 % Dirichlet boundary condition | |
| 74 case {'D','d','dirichlet'} | |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
75 tau = s * 1i*d; |
| 0 | 76 closure = obj.Hi*tau*e'; |
| 77 | |
| 78 switch class(data) | |
| 79 case 'double' | |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
80 penalty = -obj.Hi*tau*data; |
| 0 | 81 case 'function_handle' |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
82 penalty = @(t)-obj.Hi*tau*data(t); |
| 0 | 83 otherwise |
| 84 error('Wierd data argument!') | |
| 85 end | |
| 86 | |
| 87 % Unknown, boundary condition | |
| 88 otherwise | |
| 89 error('No such boundary condition: type = %s',type); | |
| 90 end | |
| 91 end | |
| 92 | |
|
946
706d1c2b4199
Raname opts to type in a bunch of interface methods
Jonatan Werpers <jonatan@werpers.com>
parents:
910
diff
changeset
|
93 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) |
| 0 | 94 % u denotes the solution in the own domain |
| 95 % v denotes the solution in the neighbour domain | |
| 96 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); | |
| 97 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 98 | |
|
67
446d67a49cd8
Fixed some errors in scheme.Schrodinger.m
Jonatan Werpers <jonatan@werpers.com>
parents:
65
diff
changeset
|
99 a = -s_u* 1/2 * 1i ; |
|
65
33f0654a2413
Fixed mistakes in Schrodinger scheme. BC are now working.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
100 b = a'; |
| 0 | 101 |
| 102 tau = b*d_u; | |
|
67
446d67a49cd8
Fixed some errors in scheme.Schrodinger.m
Jonatan Werpers <jonatan@werpers.com>
parents:
65
diff
changeset
|
103 sig = -a*e_u; |
| 0 | 104 |
| 105 closure = obj.Hi * (tau*e_u' + sig*d_u'); | |
| 106 penalty = obj.Hi * (-tau*e_v' - sig*d_v'); | |
| 107 end | |
| 108 | |
| 109 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 110 % The right boundary is considered the positive boundary | |
| 111 function [e,d,s] = get_boundary_ops(obj,boundary) | |
| 112 switch boundary | |
| 113 case 'l' | |
| 114 e = obj.e_l; | |
| 115 d = obj.d1_l; | |
| 116 s = -1; | |
| 117 case 'r' | |
| 118 e = obj.e_r; | |
| 119 d = obj.d1_r; | |
| 120 s = 1; | |
| 121 otherwise | |
| 122 error('No such boundary: boundary = %s',boundary); | |
| 123 end | |
| 124 end | |
| 125 | |
| 126 function N = size(obj) | |
| 127 N = obj.m; | |
| 128 end | |
| 129 | |
| 130 end | |
| 131 | |
| 132 methods(Static) | |
| 133 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u | |
| 134 % and bound_v of scheme schm_v. | |
| 135 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') | |
| 136 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) | |
| 137 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); | |
| 138 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); | |
| 139 end | |
| 140 end | |
| 141 end |