Mercurial > repos > public > sbplib

changeset 1021:cc61dde120cd feature/advectionRV

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add upwind dissipation to the operator inside Utux2d
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 19 Dec 2018 20:00:27 +0100
parents 5359a61cb4d9
children 234c1c02ea39
files +rv/constructDiffOps.m +scheme/Utux2d.m
diffstat 2 files changed, 22 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/constructDiffOps.m	Fri Dec 14 18:33:10 2018 +0100
+++ b/+rv/constructDiffOps.m	Wed Dec 19 20:00:27 2018 +0100
@@ -1,13 +1,14 @@
-function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs)
+function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
+    default_arg('fluxSplitting',[]);
     assert(size(BCs,1) == grid.D());
-    
+
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs);
+    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting);
     D2 = constructSymmetricD2Operator(grid, order, opSet);
     D_rv = @(v,viscosity)(D + D2(viscosity))*v;
 
     %% DiffOps for residual viscosity
-    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs);
+    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, grid, order-2, opSet, waveSpeed, BCs, fluxSplitting);
     % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
     % change the sign.
     D_flux = -D_flux;
@@ -16,12 +17,14 @@
     DvDt = @(v)D*v;
 end
 
-function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs)
-    diffOp = scheme(grid, order, opSet, waveSpeed);
+function [D, penalties] = constructTotalFluxDiffOp(scheme, grid, order, opSet, waveSpeed, BCs, fluxSplitting)
+    if isequal(opSet, @sbp.D1Upwind)
+        assert(size(fluxSplitting,1) == grid.D());
+        diffOp = scheme(grid, order, opSet, waveSpeed, fluxSplitting);
+    else
+        diffOp = scheme(grid, order, opSet, waveSpeed);
+    end
     [D, penalties]  = addClosuresToDiffOp(diffOp, BCs);
-    if (isequal(opSet,@sbp.D1Upwind))
-        D = D + getUpwindDissOp(diffOp, waveSpeed);
-    end
 end
 
 function [D, penalties] = addClosuresToDiffOp(diffOp, BCs)
@@ -35,24 +38,11 @@
     end
 end
 
-% TODO: Pass flux splitting components here instead of hard-coding the splitting.
-% TBD: Decide if this should be done in the scheme instead of by the user.
-function  DissOp = getUpwindDissOp(diffOp, waveSpeed)
-    switch diffOp.grid.D()
-        case 1
-            DissOp = -(abs(waveSpeed(1))*diffOp.Diss);
-        case 2
-            DissOp = -(abs(waveSpeed(1))*diffOp.DissOpx + abs(waveSpeed(2))*diffOp.DissOpy);
-        otherwise
-            error('3D not yet implemented')
-    end
-    
-end
-
 function D2 = constructSymmetricD2Operator(grid, order, opSet)
     % TODO: 
     % Currently only implemented for upwind operators.
     % Remove this part once the time-dependent D2 operator is implemented for other opSets
+    % or if it is decided that it should only be supported for upwind operators.
     assert(isequal(opSet,@sbp.D1Upwind))
 
     m = grid.size();
--- a/+scheme/Utux2d.m	Fri Dec 14 18:33:10 2018 +0100
+++ b/+scheme/Utux2d.m	Wed Dec 19 20:00:27 2018 +0100
@@ -27,10 +27,11 @@
 
 
     methods
-         function obj = Utux2d(g ,order, opSet, a)
+         function obj = Utux2d(g ,order, opSet, a, fluxSplitting)
 
             default_arg('a',1/sqrt(2)*[1, 1]);
             default_arg('opSet',@sbp.D2Standard);
+            default_arg('fluxSplitting',[]);
 
             assertType(g, 'grid.Cartesian');
             if iscell(a)
@@ -98,7 +99,13 @@
             obj.h = [ops_x.h ops_y.h];
             obj.order = order;
             obj.a = a;
-            obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
+
+            if (isequal(opSet,@sbp.D1Upwind))
+                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy + fluxSplitting{1}*obj.DissOpx + fluxSplitting{2}*obj.DissOpy);
+            else
+                obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
+            end
+
         end
         % Closure functions return the opertors applied to the own domain to close the boundary
         % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.