diff +noname/calculateErrors.m @ 832:5573913a0949 feature/burgers1d

Merged with default, and updated +scheme/Burgers1D accordingly
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 11 Sep 2018 15:58:35 +0200
parents 1201eb16557e
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+noname/calculateErrors.m	Tue Sep 11 15:58:35 2018 +0200
@@ -0,0 +1,53 @@
+% [discr, trueSolution] =  schemeFactory(m)
+%     where trueSolution should be a timeSnapshot of the true solution a time T
+% T is the end time
+% m are grid size parameters.
+% N are number of timesteps to use for each gird size
+% timeOpt are options for the timeStepper
+% errorFun is a function_handle taking 2 or 3 arguments, errorFun(trueSolution, approxSolution), errorFun(trueSolution, approxSolution, discr)
+function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt)
+    %TODO: Ability to choose paralell or not
+    assertType(schemeFactory, 'function_handle');
+    assertNumberOfArguments(schemeFactory, 1);
+    assertScalar(T);
+    assert(length(m) == length(N), 'Vectors m and N must have the same length');
+    assertType(errorFun, 'function_handle');
+
+    if ~ismember(nargin(errorFun), [2,3])
+        error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun));
+    end
+
+    default_arg('timeOpt', struct());
+
+
+    e = zeros(1,length(m));
+    parfor i = 1:length(m)
+        done = timeTask('m = %3d ', m(i));
+
+        [discr, trueSolution] = schemeFactory(m(i));
+
+        timeOptTemp = timeOpt;
+        timeOptTemp.k = T/N(i);
+        ts = discr.getTimestepper(timeOptTemp);
+        ts.stepTo(N(i), true);
+        approxSolution = discr.getTimeSnapshot(ts);
+
+        switch nargin(errorFun)
+            case 2
+                e(i) = errorFun(trueSolution, approxSolution);
+            case 3
+                e(i) = errorFun(trueSolution, approxSolution, discr);
+        end
+
+        fprintf('e = %.4e', e(i))
+        done()
+    end
+    fprintf('\n')
+end
+
+
+%% Example error function
+% u_true = grid.evalOn(dr.grid, @(x,y)trueSolution(T,x,y));
+% err = u_true-u_false;
+% e(i) = norm(err)/norm(u_true);
+% % e(i) = sqrt(err'*d.H*d.J*err/(u_true'*d.H*d.J*u_true));