Mercurial > repos > public > sbplib

changeset 1156:8c0e2b50f018 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt to fix oscilations when not using bdfs by adding lower-order 2nd derivative correction term to the Residual
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 25 Jun 2019 13:06:42 +0200
parents 336ee37a0617
children 9ad15100270d
files +rv/constructDiffOps.m
diffstat 1 files changed, 35 insertions(+), 9 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/constructDiffOps.m	Mon Jun 24 17:44:21 2019 +0200
+++ b/+rv/constructDiffOps.m	Tue Jun 25 13:06:42 2019 +0200
@@ -1,16 +1,43 @@
-function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs)
+function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, schemeOrder, residualOrder, schemeParams, opSet, BCs)
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
-    D2 = constructSymmetricD2Operator(g, order, opSet);
+    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs);
+    D2 = constructSymmetricD2Operator(g, schemeOrder, opSet);
     D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v);
 
     %% DiffOps for residual viscosity
-    [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
-    % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
+    [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, residualOrder, schemeParams, opSet, BCs);
+    % TODO: Construct D_flux without closures when using bdfs.
+    % diffOp = scheme(g, residualOrder, schemeParams{:}, opSet);
+    % if ~isa(diffOp.D, 'function_handle')
+    %     D_flux = @(v) diffOp.D*v;
+    % else
+    %     D_flux = diffOp.D;
+    % end
+
+    % DiffOp for flux in residual viscosity. Due to sign conventions of the implemented schemes, we need to
     % change the sign.
     D_flux = @(v) -D_flux(v);
+    lim = {g.x{1}(1), g.x{1}(end)};
+    o = 2;
+    h = g.scaling();
+    gam = 1;
+    if isequal(opSet,@sbp.D1Upwind)
+        ops = sbp.D1Upwind(g.m(), lim, o);
+        Hi = ops.HI;
+        e_r = ops.e_r;
+        e_l = ops.e_l;
+        B = e_r*e_r' - e_l*e_l';
+        M =  ops.Dm - Hi*B;
+        D2 = M*ops.Dp;
+    else
+        ops = sbp.D2Standard(g.m(), lim, o);
+        D2 = ops.D2;
+    end
+    D2_low = gam*h^(residualOrder)*D2;
+    D2_low = @(v) D2_low*v;
     % DiffOp for time derivative in residual viscosity
-    DvDt = D;
+    DvDt = @(v) D(v)+D2_low(v);
+    %DvDt = @(v) D(v);
 end
 
 function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
@@ -37,8 +64,6 @@
 end
 
 function D2 = constructSymmetricD2Operator(g, order, opSet)
-
-
     m = g.size();
     ops = cell(g.D(),1);
     I = cell(g.D(),1);
@@ -61,7 +86,8 @@
                 D2 = @(Viscosity) M*Viscosity*Dp;
             else
                 % TODO: Fix closure for D2Variable
-                D2 = @(Viscosity)ops{1}.D2(Viscosity);
+                % TODO: Fix Viscosity not being vector
+                D2 = @(Viscosity)ops{1}.D2(diag(Viscosity));
             end
         case 2
             % TODO: