Mercurial > repos > public > sbplib

comparison +rv/+diffops/constructDiffOps.m @ 1191:881afc40a3d2 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change to use the interpolation operators in sbplib
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 29 Jul 2019 16:43:08 +0200
parents 79618b58b0a0
children 5271c4670733
comparison
equal deleted inserted replaced
1190:79618b58b0a0 1191:881afc40a3d2
6 diffOpsStruct = diffOpsBdf(varargin{:}); 6 diffOpsStruct = diffOpsBdf(varargin{:});
7 case {'ms','multi-stage'} 7 case {'ms','multi-stage'}
8 diffOpsStruct = diffOpsMultiStage(varargin{:}); 8 diffOpsStruct = diffOpsMultiStage(varargin{:});
9 case {'mg1','multi-grid1'} 9 case {'mg1','multi-grid1'}
10 diffOpsStruct = diffOpsMultiGrid1(varargin{:}); 10 diffOpsStruct = diffOpsMultiGrid1(varargin{:});
11 case {'mg2','multi-grid2'} 11 case {'mg','mg2','multi-grid2'} % Default multigrid diffops
12 diffOpsStruct = diffOpsMultiGrid2(varargin{:}); 12 diffOpsStruct = diffOpsMultiGrid2(varargin{:});
13 end 13 end
14 end 14 end
15 15
16 function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) 16 function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
83 diffOpsStruct = struct('D_scheme', D_scheme,... 83 diffOpsStruct = struct('D_scheme', D_scheme,...
84 'D_flux', D_flux,... 84 'D_flux', D_flux,...
85 'D_t', D_t,... 85 'D_t', D_t,...
86 'penalties_scheme', {penalties_scheme}); 86 'penalties_scheme', {penalties_scheme});
87 end 87 end
88 %% Multigrid functions:
89 % TODO: Implement properly.
90 function v_c = coarse_v(g, g_c, v)
91 V = grid.funcToMatrix(g,v);
92 v_c = reshape(V(1:2:end,1:2:end)',g_c.N(),1);
93 end
94 88
95 % Note: Only works for equidistant grids. 89 % Note: Only works for equidistant grids.
96 % TODO: Change from using matlabs interp to using proper interpolation operators..
97 function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) 90 function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
98 m = g.m(); 91 m = g.m();
99 lim = g.lim(); 92 lim = g.lim();
100 m_c = (m-1)/2 + 1; 93 m_c = (m-1)/2 + 1;
101 switch g.D() 94 switch g.D()
102 case 1 95 case 1
103 g_c = grid.equidistant(m_c, lim{1}); 96 interpOps = sbp.InterpOpsOP(m_c, m, schemeOrder, schemeOrder);
104 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); 97 I = interpOps.Iu2v.good;
105 x = g.x{1}; 98 g_c = grid.equidistant(m_c, lim{1});
106 x_c = x(1:2:end); 99 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
107 D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline'); 100 D_c = @(v) I*D_c(v(1:2:end));
108 case 2 101 case 2
109 g_c = grid.equidistant(m_c, lim{1}, lim{2}); 102 interpOps_x = sbp.InterpOpsOP(m_c(1), m(1), schemeOrder, schemeOrder);
110 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); 103 I_x = interpOps_x.Iu2v.good;
111 D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1); 104 interpOps_y = sbp.InterpOpsOP(m_c(2), m(2), schemeOrder, schemeOrder);
105 I_y = interpOps_y.Iu2v.good;
106 I = kron(I_x,I_y);
107 g_c = grid.equidistant(m_c, lim{1}, lim{2});
108 ind = grid.funcToMatrix(g, 1:g.N());
109 ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
110 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
111 D_c = @(v) I*D_c(v(ind_c));
112 end 112 end
113 end 113 end
114 114
115 function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs) 115 function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs)
116 m = g.m(); 116 m = g.m();
124 x = g.x{1}; 124 x = g.x{1};
125 x_c = x(1:2:end); 125 x_c = x(1:2:end);
126 D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline'); 126 D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
127 case 2 127 case 2
128 g_c = grid.equidistant(m_c, lim{1}, lim{2}); 128 g_c = grid.equidistant(m_c, lim{1}, lim{2});
129 ind = grid.funcToMatrix(g, 1:g.N());
130 ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
129 [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs); 131 [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
130 D_c = rv.diffops.addClosuresToDiffOp(D_c, closures); 132 D_c = rv.diffops.addClosuresToDiffOp(D_c, closures);
131 D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1); 133 D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(v(ind_c))),'spline')',g.N(),1);
132 end 134 end
133 end 135 end
134 136