Mercurial > repos > public > sbplib

annotate +multiblock/stitchSchemes.m @ 692:dd3e9a0f5032 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added jacobian to stichshemes
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 15 Sep 2017 09:21:29 +0200
parents 324c927d8b1d
children a5d4e3ced9a5
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 % Stitch schemes together given connection matrix and BC vector.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 % schmHand - function_handle to a Scheme constructor
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 % order - order of accuracy
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 % schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 % blocks - block definitions, On whatever form the scheme expects.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 % ms - grid points in each direction for each block. Ex {[10,10], [10, 20]}
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 % conn - connection matrix
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 % bound - boundary condition vector, array of structs with fields w,e,s,n
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 % each field with a parameter array that is sent to schm.boundary_condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 %
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 % Output parameters are cell arrays and cell matrices.
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
12 %
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
13 % Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound)
692
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
14 function [schms, D, H,penalty,J] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep)
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
15 default_arg('schmParam',[]);
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
16 default_arg('timeDep','N');
181
419ec303e97d Made some local changes to stichSchemes.
Jonatan Werpers <jonatan@werpers.com>
parents: 120
diff changeset
17 n_blocks = numel(grids);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 % Creating Schemes
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 for i = 1:n_blocks
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
21 if isempty(schmParam);
181
419ec303e97d Made some local changes to stichSchemes.
Jonatan Werpers <jonatan@werpers.com>
parents: 120
diff changeset
22 schms{i} = schmHand(grids{i},order,[]);
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
23 elseif ~iscell(schmParam)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 param = schmParam(i);
181
419ec303e97d Made some local changes to stichSchemes.
Jonatan Werpers <jonatan@werpers.com>
parents: 120
diff changeset
25 schms{i} = schmHand(grids{i},order,param);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 else
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 param = schmParam{i};
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
28 if iscell(param)
181
419ec303e97d Made some local changes to stichSchemes.
Jonatan Werpers <jonatan@werpers.com>
parents: 120
diff changeset
29 schms{i} = schmHand(grids{i},order,param{:});
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
30 else
181
419ec303e97d Made some local changes to stichSchemes.
Jonatan Werpers <jonatan@werpers.com>
parents: 120
diff changeset
31 schms{i} = schmHand(grids{i},order,param);
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
32 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 % class(schmParam)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 % class(ms)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 % class(blocks)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 % class(schmParam{i})
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 % class(ms)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40
76
5c569cbef49e Added more input parameter handling and fixed some bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
41
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 % Total norm
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 H = cell(n_blocks,n_blocks);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 for i = 1:n_blocks
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 H{i,i} = schms{i}.H;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 end
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
50
692
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
51 J = cell(n_blocks,n_blocks);
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
52 for i = 1:n_blocks
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
53 J{i,i} = schms{i}.J;
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
54 end
dd3e9a0f5032 added jacobian to stichshemes
Ylva Rydin <ylva.rydin@telia.com>
parents: 498
diff changeset
55
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 %% Total system matrix
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
57
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 % Differentiation terms
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 D = cell(n_blocks,n_blocks);
498
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
60 penalty = cell(n_blocks,n_blocks);
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
61
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
62
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 for i = 1:n_blocks
498
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
64 penalty{i,i}= @(t)0;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 D{i,i} = schms{i}.D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 end
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
67
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 % Boundary penalty terms
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
69 switch timeDep
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
70 case {'n','no','N','No'}
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
71 for i = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
72 if ~isstruct(bound{i})
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
73 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
74 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
75
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
76 fn = fieldnames(bound{i});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
77 for j = 1:length(fn)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
78 bc = bound{i}.(fn{j});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
79 if isempty(bc)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
80 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
81 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
82
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
83 [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
84 D{i,i} = D{i,i}+closure;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
85 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 end
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
87
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
88 % Interface penalty terms
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
89 for i = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
90 for j = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
91 intf = conn{i,j};
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
92 if isempty(intf)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
93 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
94 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
95
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
96 [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
97 D{i,i} = D{i,i} + uu;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
98 D{i,j} = uv;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
99 D{j,j} = D{j,j} + vv;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
100 D{j,i} = vu;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
101 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 end
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
103 case {'y','yes','Y','Yes'}
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
104 for i = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
105 if ~isstruct(bound{i})
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
106 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
107 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
108
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
109 fn = fieldnames(bound{i});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
110 for j = 1:length(fn)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
111 bc = bound{i}.(fn{j});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
112 if isempty(bc)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
113 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
114 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
115
498
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
116 [closure, penaltyi] = schms{i}.boundary_condition(fn{j},bc{:});
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
117 D{i,i} =@(t) D{i,i}(t) + closure(t);
498
324c927d8b1d chnaged sbp interfacein 1d among many things
Ylva Rydin <ylva.rydin@telia.com>
parents: 496
diff changeset
118 penalty{i,i} = @(t) penalty{i,i}(t) + penaltyi(t);
496
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
119 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
120 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
121
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
122 % Interface penalty terms
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
123 for i = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
124 for j = 1:n_blocks
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
125 intf = conn{i,j};
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
126 if isempty(intf)
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
127 continue
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
128 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
129
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
130 [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
131 D{i,i} = @(t) D{i,i}(t) + uu(t);
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
132 D{i,j} = uv;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
133 D{j,j} = @(t)D{j,j}(t) + vv(t);
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
134 D{j,i} = vu;
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
135 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
136 end
437fad4a47e1 Add 1d interface for shrodinger in 1d
Ylva Rydin <ylva.rydin@telia.com>
parents: 181
diff changeset
137 end