Mercurial > repos > public > sbplib

comparison +rv/+diffops/constructDiffOps.m @ 1184:ecc605453733 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move specialized constructDiffOps* into constructDiffOps.m. Add multigrid RV using the standard RungekuttaRV time stepper
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 05 Jul 2019 17:51:11 +0200
parents e3d8f24b2c1c
children 79618b58b0a0
comparison
equal deleted inserted replaced
1183:27f3b173dcfa 1184:ecc605453733
1 function diffOpsStruct = constructDiffOps(method, varargin) 1 function diffOpsStruct = constructDiffOps(method, varargin)
2 switch method 2 switch method
3 case {'bdf','backwards-difference-formula'} 3 case {'standard'}
4 diffOpsStruct = rv.diffops.constructDiffOpsBdf(varargin{:}); 4 diffOpsStruct = diffOpsStandard(varargin{:});
5 case {'ms','multi-stage'} 5 case {'bdf','backwards-difference-formula'}
6 diffOpsStruct = rv.diffops.constructDiffOpsMultiStage(varargin{:}); 6 diffOpsStruct = diffOpsBdf(varargin{:});
7 case {'mg','multi-grid'} 7 case {'ms','multi-stage'}
8 diffOpsStruct = rv.diffops.constructDiffOpsMultiGrid(varargin{:}); 8 diffOpsStruct = diffOpsMultiStage(varargin{:});
9 case {'mg1','multi-grid1'}
10 diffOpsStruct = diffOpsMultiGrid1(varargin{:});
11 case {'mg2','multi-grid2'}
12 diffOpsStruct = diffOpsMultiGrid2(varargin{:});
13 end
9 end 14 end
15
16 function diffOpsStruct = diffOpsStandard(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
17 % DiffOps for stabilized solution vector
18 [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
19 %% DiffOps for residual viscosity
20 D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
21 D_flux = @(v) -D_flux(v);
22 diffOpsStruct = struct('D_scheme', D_scheme,...
23 'D_flux', D_flux,...
24 'D_t', D_t,...
25 'penalties_scheme', penalties_scheme);
26 end
27
28 function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
29 %% DiffOps for solution vector
30 [D_scheme, penalties_scheme, ~] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
31 %% DiffOps for residual viscosity
32 [D_flux, ~] = rv.diffops.constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
33 D_flux = @(v) -D_flux(v);
34 diffOpsStruct = struct('D_scheme', D_scheme,...
35 'D_flux', D_flux,...
36 'penalties_scheme', penalties_scheme);
37 end
38
39 function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
40 % DiffOps for stabilized solution vector
41 [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
42 % DiffOp for unstabilized solution vector
43 D_unstable = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs);
44 %% DiffOps for residual viscosity
45 [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
46 D_flux = @(v) -D_t(v);
47 diffOpsStruct = struct('D_scheme', D_scheme,...
48 'D_unstable', D_unstable,...
49 'D_flux', D_flux,...
50 'D_t', D_t,...
51 'penalties_scheme', penalties_scheme,...
52 'penalties_res', penalties_res);
53 end
54
55 function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
56 % DiffOps for stabilized solution vector
57 [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
58 % DiffOp for unstabilized solution vector
59 D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
60 %% DiffOps for residual viscosity
61 [D_flux, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
62 D_flux = @(v) -D_flux(v);
63 D_t = @(v) D(v);
64 diffOpsStruct = struct('D_scheme', D_scheme,...
65 'D_coarse', D_coarse,...
66 'D_flux', D_flux,...
67 'D_t', D_t,...
68 'penalties_scheme', penalties_scheme,...
69 'penalties_res', penalties_res);
70 end
71
72 function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
73 % DiffOps for stabilized solution vector
74 [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
75 % TODO: What orders to use here?
76 D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
77 %% DiffOps for residual viscosity
78 D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs);
79 D_flux = @(v) -D_flux(v);
80 D_t = @(v) D_coarse(v);
81 diffOpsStruct = struct('D_scheme', D_scheme,...
82 'D_flux', D_flux,...
83 'D_t', D_t,...
84 'penalties_scheme', penalties_scheme);
85 end
86
87 % TODO: Only works for equidistant grids
88 function D_coarse = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
89 m = g.m();
90 lim = g.lim();
91 m_coarse = (m-1)/2 + 1;
92 g_coarse = grid.equidistant(m_coarse, lim{1});
93 D_coarse = rv.diffops.constructFluxDiffOpWithClosures(scheme, g_coarse, schemeOrder, schemeParams, opSet, BCs);
94 x = g.x{1};
95 x_coarse = x(1:2:end);
96 % TODO: Fix interpolation
97 D_coarse = @(v) interp1(x_coarse,D_coarse(v(1:2:end)),x,'spline');
98 end