Mercurial > repos > public > sbplib

changeset 1151:03ecf18d035f feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Set RV artificial viscosity to zero around borders of the domain
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 18 Feb 2019 09:00:00 +0100
parents 922996695952
children 010bb2677230
files +rv/ResidualViscosity.m
diffstat 1 files changed, 17 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/ResidualViscosity.m	Thu Jan 24 10:23:52 2019 +0100
+++ b/+rv/ResidualViscosity.m	Mon Feb 18 09:00:00 2019 +0100
@@ -7,20 +7,21 @@
         h % Length scale used for scaling the viscosity. Typically grid spacing.
         normalization % Function used to normalize the residual such that it is amplified in the
                       % shocks.
+        boundaryIndices % indices of the boundary of the domain
     end
 
     methods
         %  TBD: Decide on how to treat waveSpeed. It would be nice to just pass a constant value without
         %       wrapping it in a function.
-        function obj = ResidualViscosity(Df, waveSpeed, Cmax, Cres, h, normalization)
+        function obj = ResidualViscosity(g, Df, waveSpeed, Cmax, Cres, h, normalization)
             default_arg('normalization',@(v)abs(obj.minmaxDiffNeighborhood1d(v)-norm(v-mean(v),inf)));
-            %default_arg('normalization',@(v)norm(v/2,inf));
             obj.Df = Df;
             obj.waveSpeed = waveSpeed;
             obj.h = h;
             obj.Cmax = Cmax;
             obj.Cres = Cres;
             obj.normalization = normalization;
+            obj.boundaryIndices = obj.getBoundaryIndices(g);
         end
 
         function [viscosity, Df, firstOrderViscosity, residualViscosity] = evaluate(obj, v, dvdt)
@@ -28,10 +29,24 @@
             firstOrderViscosity = obj.Cmax*obj.h*abs(obj.waveSpeed(v));
             residualViscosity = obj.Cres*obj.h^2*abs(dvdt + Df)./obj.normalization(v);
             viscosity = min(firstOrderViscosity, residualViscosity);
+            % TODO: If the boundary conditions are implemented correctly, this might not be needed.
+            viscosity(obj.boundaryIndices) = 0;
         end
     end
     methods (Static)
 
+        function boundaryind = getBoundaryIndices(g)
+            ind = grid.funcToMatrix(g, 1:g.N()); 
+            switch g.D()
+            case 1
+                boundaryind = [ind(1) ind(end)];
+            case 2
+                boundaryind = [ind(1,:) ind(end,:), ind(:,1)' ind(:,end)'];
+            case 3
+                error;
+            end
+        end
+
         % function f = minmaxDiffNeighborhood(g)
         %     switch g.D()
         %         case 1