Mercurial > repos > public > sbplib

comparison +rv/+diffops/constructDiffOps.m @ 1190:79618b58b0a0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor constructDiffops and remove closures from some of the diffops used to compute the residual.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 08 Jul 2019 15:12:19 +0200
parents ecc605453733
children 881afc40a3d2
comparison
equal deleted inserted replaced
1189:6cb03209f0a7 1190:79618b58b0a0
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)
17 % DiffOps for stabilized solution vector 17 % DiffOps for stabilized solution vector
18 [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); 18 [D_scheme, penalties_scheme, D_t] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
19 %% DiffOps for residual viscosity 19 %% DiffOps for residual viscosity
20 D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); 20 D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
21 D_flux = @(v) -D_flux(v); 21 D_flux = @(v) -D_flux(v);
22 diffOpsStruct = struct('D_scheme', D_scheme,... 22 diffOpsStruct = struct('D_scheme', D_scheme,...
23 'D_flux', D_flux,... 23 'D_flux', D_flux,...
24 'D_t', D_t,... 24 'D_t', D_t,...
25 'penalties_scheme', penalties_scheme); 25 'penalties_scheme', {penalties_scheme});
26 end 26 end
27 27
28 function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) 28 function diffOpsStruct = diffOpsBdf(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
29 %% DiffOps for solution vector 29 %% DiffOps for solution vector
30 [D_scheme, penalties_scheme, ~] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); 30 [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
31 %% DiffOps for residual viscosity 31 %% DiffOps for residual viscosity
32 [D_flux, ~] = rv.diffops.constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); 32 D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
33 D_flux = @(v) -D_flux(v); 33 D_flux = @(v) -D_flux(v);
34 diffOpsStruct = struct('D_scheme', D_scheme,... 34 diffOpsStruct = struct('D_scheme', D_scheme,...
35 'D_flux', D_flux,... 35 'D_flux', D_flux,...
36 'penalties_scheme', penalties_scheme); 36 'penalties_scheme', {penalties_scheme});
37 end 37 end
38 38
39 % TODO: Remove?
39 function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) 40 function diffOpsStruct = diffOpsMultiStage(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
40 % DiffOps for stabilized solution vector 41 % DiffOps for stabilized solution vector
41 [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); 42 [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
42 % DiffOp for unstabilized solution vector 43 % DiffOp for unstabilized solution vector
43 D_unstable = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); 44 [D_unstable, closures] = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
45 D_unstable = rv.diffops.addClosuresToDiffOp(D_unstable, closures);
44 %% DiffOps for residual viscosity 46 %% DiffOps for residual viscosity
45 [D_t, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); 47 D_t = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
46 D_flux = @(v) -D_t(v); 48 D_flux = @(v) -D_t(v);
49 % TODO: Use residual penalties as well here?
47 diffOpsStruct = struct('D_scheme', D_scheme,... 50 diffOpsStruct = struct('D_scheme', D_scheme,...
48 'D_unstable', D_unstable,... 51 'D_unstable', D_unstable,...
49 'D_flux', D_flux,... 52 'D_flux', D_flux,...
50 'D_t', D_t,... 53 'D_t', D_t,...
51 'penalties_scheme', penalties_scheme,... 54 'penalties_scheme', {penalties_scheme});
52 'penalties_res', penalties_res);
53 end 55 end
54 56
55 function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) 57 function diffOpsStruct = diffOpsMultiGrid1(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
56 % DiffOps for stabilized solution vector 58 % DiffOps for stabilized solution vector
57 [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); 59 [D_scheme, penalties_scheme, D_f] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
58 % DiffOp for unstabilized solution vector 60 % DiffOp for unstabilized solution vector
59 D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); 61 D_coarse = coarserGridDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs);
60 %% DiffOps for residual viscosity 62 %% DiffOps for residual viscosity
61 [D_flux, penalties_res] = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, residualOrder, schemeParams, opSet, BCs); 63 D_flux = rv.diffops.constructFluxDiffOps(scheme, g, residualOrder, schemeParams, opSet, BCs);
62 D_flux = @(v) -D_flux(v); 64 D_flux = @(v) -D_flux(v);
63 D_t = @(v) D(v); 65 D_t = @(v) D_f(v);
66 % TODO: Use residual penalties as well here?
64 diffOpsStruct = struct('D_scheme', D_scheme,... 67 diffOpsStruct = struct('D_scheme', D_scheme,...
65 'D_coarse', D_coarse,... 68 'D_coarse', D_coarse,...
66 'D_flux', D_flux,... 69 'D_flux', D_flux,...
67 'D_t', D_t,... 70 'D_t', D_t,...
68 'penalties_scheme', penalties_scheme,... 71 'penalties_scheme', {penalties_scheme});
69 'penalties_res', penalties_res);
70 end 72 end
71 73
72 function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs) 74 function diffOpsStruct = diffOpsMultiGrid2(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
73 % DiffOps for stabilized solution vector 75 % DiffOps for stabilized solution vector
74 [D_scheme, penalties_scheme, D] = rv.diffops.constructSchemeDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs); 76 [D_scheme, penalties_scheme] = rv.diffops.constructSchemeDiffOps(scheme, g, schemeOrder, schemeParams, opSet, BCs);
75 % TODO: What orders to use here? 77 % TODO: What orders to use here?
76 D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs); 78 D_coarse = coarserGridDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
77 %% DiffOps for residual viscosity 79 %% DiffOps for residual viscosity
78 D_flux = rv.diffops.constructFluxDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs); 80 D_flux = rv.diffops.constructFluxDiffOps(scheme, g, schemeOrder , schemeParams, opSet, BCs);
79 D_flux = @(v) -D_flux(v); 81 D_flux = @(v) -D_flux(v);
80 D_t = @(v) D_coarse(v); 82 D_t = @(v) D_coarse(v);
81 diffOpsStruct = struct('D_scheme', D_scheme,... 83 diffOpsStruct = struct('D_scheme', D_scheme,...
82 'D_flux', D_flux,... 84 'D_flux', D_flux,...
83 'D_t', D_t,... 85 'D_t', D_t,...
84 'penalties_scheme', penalties_scheme); 86 'penalties_scheme', {penalties_scheme});
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);
85 end 93 end
86 94
87 % TODO: Only works for equidistant grids 95 % Note: Only works for equidistant grids.
88 function D_coarse = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs) 96 % TODO: Change from using matlabs interp to using proper interpolation operators..
97 function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
89 m = g.m(); 98 m = g.m();
90 lim = g.lim(); 99 lim = g.lim();
91 m_coarse = (m-1)/2 + 1; 100 m_c = (m-1)/2 + 1;
92 g_coarse = grid.equidistant(m_coarse, lim{1}); 101 switch g.D()
93 D_coarse = rv.diffops.constructFluxDiffOpWithClosures(scheme, g_coarse, schemeOrder, schemeParams, opSet, BCs); 102 case 1
94 x = g.x{1}; 103 g_c = grid.equidistant(m_c, lim{1});
95 x_coarse = x(1:2:end); 104 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
96 % TODO: Fix interpolation 105 x = g.x{1};
97 D_coarse = @(v) interp1(x_coarse,D_coarse(v(1:2:end)),x,'spline'); 106 x_c = x(1:2:end);
107 D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
108 case 2
109 g_c = grid.equidistant(m_c, lim{1}, lim{2});
110 D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
111 D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1);
112 end
98 end 113 end
114
115 function D_c = coarserGridDiffOpWithClosures(scheme, g, schemeOrder, schemeParams, opSet, BCs)
116 m = g.m();
117 lim = g.lim();
118 m_c = (m-1)/2 + 1;
119 switch g.D()
120 case 1
121 g_c = grid.equidistant(m_c, lim{1});
122 [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
123 D_c = rv.diffops.addClosuresToDiffOp(D_c, closures);
124 x = g.x{1};
125 x_c = x(1:2:end);
126 D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
127 case 2
128 g_c = grid.equidistant(m_c, lim{1}, lim{2});
129 [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
130 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);
132 end
133 end
134