Mercurial > repos > public > sbplib
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));