From 6c7ae1ce1878b11e0b49bec06f642aae4a4aa090 Mon Sep 17 00:00:00 2001 From: David Green Date: Tue, 26 May 2020 10:33:52 -0400 Subject: [PATCH 1/4] more work towards BDF2 --- runtime_defaults.m | 4 +-- time_advance.m | 64 ++++++++++++++++++++++++++++++---- two_scale/two_scale_rel_2.mat | Bin 413 -> 413 bytes 3 files changed, 59 insertions(+), 9 deletions(-) diff --git a/runtime_defaults.m b/runtime_defaults.m index 022a8ac5..b992809a 100644 --- a/runtime_defaults.m +++ b/runtime_defaults.m @@ -16,7 +16,7 @@ default_use_oldhash = false; default_use_oldcoeffmat = false; default_timestep_method = 'RK3'; -valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode45','RK3'}; +valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode45','RK3','BDF2'}; check_timestep_method = @(x) any(validatestring(x,valid_timestep_methods)); default_time_independent_A = false; default_many_solution_capable = false; @@ -125,7 +125,7 @@ opts.timestep_method = input_parser.Results.timestep_method; opts.build_A = false; -if sum(strcmp(opts.timestep_method,{'BE','CN'}))>0 +if sum(strcmp(opts.timestep_method,{'BE','CN','BDF2'}))>0 opts.build_A = true; end diff --git a/time_advance.m b/time_advance.m index 07a9c15c..0eccc535 100644 --- a/time_advance.m +++ b/time_advance.m @@ -1,17 +1,20 @@ -function f = time_advance(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax) +function f1 = time_advance(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax) if strcmp(opts.timestep_method,'BE') % Backward Euler (BE) first order - f = backward_euler(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax); + f1 = backward_euler(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); elseif strcmp(opts.timestep_method,'CN') % Crank Nicolson (CN) second order - f = crank_nicolson(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax); + f1 = crank_nicolson(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); +elseif strcmp(opts.timestep_method,'BDF2') + % Backward Difference Forumla (BDF) second order + f1 = bdf2(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); elseif sum(strcmp(opts.timestep_method,{'ode15i','ode15s','ode45'}))>0 % One of the atlab ODE integrators. - f = ODEm(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax); + f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); else % RK3 TVD - f = RungeKutta3(pde,opts,A_data,f,t,dt,deg,hash_table,Vmax,Emax); + f1 = RungeKutta3(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); end end @@ -49,7 +52,7 @@ if strcmp(opts.timestep_method,'ode45') if(~opts.quiet);disp('Using ode45');end - options = odeset('RelTol',1e-3,'AbsTol',1e-6,'Stats','off'); + options = odeset('RelTol',1e-3,'AbsTol',1e-6,'Stats','on'); [tout,fout] = ode45(@explicit_ode,[t0 t0+dt],f0,options); elseif strcmp(opts.timestep_method,'ode15s') @@ -72,7 +75,7 @@ % call ode15s if(~opts.quiet);disp('Using ode15s');end - options = odeset('RelTol',1e-4,'AbsTol',1e-3,... + options = odeset('RelTol',1e-4,'AbsTol',1e-6,... 'Stats','on','OutputFcn',@odetpbar,'Refine',20);%,'Jacobian', J2);%'JPattern',S); [tout,fout] = ode15s(@explicit_ode,[t0 t0+dt],f0,options); @@ -237,3 +240,50 @@ end end + + +function f2 = bdf2(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax) + +s1 = source_vector(pde,opts,hash_table,t+dt); +s2 = source_vector(pde,opts,hash_table,t+2*dt); + +bc1 = boundary_condition_vector(pde,opts,hash_table,t+dt); +bc2 = boundary_condition_vector(pde,opts,hash_table,t+2*dt); + +% %% +% Apply any non-identity LHS mass matrix coefficient + +applyLHS = ~isempty(pde.termsLHS); + +%Obtain f1 through same process as Backward Euler +if t == 0 + %opts.timestep_method = 'ode15s'; + %f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); + fac = 100; + for n=1:fac + f0 = backward_euler(pde,opts,A_data,f0,t,dt/fac,deg,hash_table,Vmax,Emax); + t = t + dt/fac; + end + f1 = f0; + save('f1.mat', 'f1'); +else + load('f1.mat', 'f1'); %Need to update this at the end +end + +if applyLHS %Now we only need BDF2 + %Obtain updated A and ALHS matrices + [~,A2,ALHS2] = apply_A(pde,opts,A_data,f1,deg); + I = eye(numel(diag(A2))); + AA2 = I - 2/3*dt*(ALHS2 \ A2); + b2 = 4*f1/3 - f0/3 + 2/3*dt*(ALHS2 \ (bc2 + s2)); +else + [~,A2] = apply_A(pde,opts,A_data,f1,deg); + I = eye(numel(diag(A2))); + AA2 = I -2/3*dt*A2; + b2 = 4*f1/3 - f0/3 + 2/3*dt*(bc2 + s2); +end +f2 = AA2 \ b2; +f1 = f2; %Updated initial guess for f1 +save('f1', 'f1'); + +end \ No newline at end of file diff --git a/two_scale/two_scale_rel_2.mat b/two_scale/two_scale_rel_2.mat index be26675300ee3321d8c4587d9f7f84dd0031c344..9422cf561b2ca4c427eebc949bbfb06ea69d0513 100644 GIT binary patch delta 61 zcmbQsJePTbypFGCUxcHXiH?GEQEFmIYKlUBo|S@cex8C$YO;c{fr6p2m5GIwv9W@Y RfuZHZK;?-EtQ$)-838V45j6k+ From 397d48310782d3a45addc5c01f9f4b1781bb3325 Mon Sep 17 00:00:00 2001 From: lwonnell Date: Tue, 26 May 2020 19:44:57 -0400 Subject: [PATCH 2/4] Added Runge-Kutta option for t==0 in time_advance and added test to asgard_test.m for BDF2 file --- asgard_test.m | 14 ++++++++++++++ time_advance.m | 15 ++++++++------- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/asgard_test.m b/asgard_test.m index 16180898..f9649baf 100644 --- a/asgard_test.m +++ b/asgard_test.m @@ -265,3 +265,17 @@ function asgard_fokkerplanck1_5p1a_implicit_test(testCase) verifyLessThan(testCase,err,2.5e-2); end +function asgard_BDF2_test(testCase) +addpath(genpath(pwd)); +disp('Testing BDF2'); +[err,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... + 'lev',3,'num_steps',2,'dt',0.001,'quiet',true); +[err1,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... + 'lev',3,'num_steps',4,'dt',0.0005,'quiet',true); +[err2,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... + 'lev',3,'num_steps',8,'dt',0.00025,'quiet',true); +logslope = (log(err1) - log(err))/(log(1/0.001) - log(1/0.0005)); +slopeerr = abs(logslope - 2.000); +verifyLessThan(testCase,slopeerr,1e-3); +end + diff --git a/time_advance.m b/time_advance.m index 0eccc535..f49bba29 100644 --- a/time_advance.m +++ b/time_advance.m @@ -256,13 +256,14 @@ applyLHS = ~isempty(pde.termsLHS); %Obtain f1 through same process as Backward Euler -if t == 0 - %opts.timestep_method = 'ode15s'; - %f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); - fac = 100; +if t <= dt +% opts.timestep_method = 'ode45'; +% f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); + fac = 200; for n=1:fac f0 = backward_euler(pde,opts,A_data,f0,t,dt/fac,deg,hash_table,Vmax,Emax); - t = t + dt/fac; + %f0 = RungeKutta3(pde,opts,A_data,f0,t,dt/fac,deg,hash_table,Vmax,Emax); + t = t + dt/fac; end f1 = f0; save('f1.mat', 'f1'); @@ -283,7 +284,7 @@ b2 = 4*f1/3 - f0/3 + 2/3*dt*(bc2 + s2); end f2 = AA2 \ b2; -f1 = f2; %Updated initial guess for f1 -save('f1', 'f1'); +f1 = f2 ; +save('f1.mat', 'f1'); end \ No newline at end of file From 004aeac382a749c4fdabb7de3aa77d64e267e63a Mon Sep 17 00:00:00 2001 From: lwonnell Date: Tue, 26 May 2020 20:12:17 -0400 Subject: [PATCH 3/4] Modified asgard_test to include degree. Diffusion1 now gives close to 2nd order accuracy for BDF2 but not advection1 --- asgard_test.m | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/asgard_test.m b/asgard_test.m index f9649baf..cd982d2e 100644 --- a/asgard_test.m +++ b/asgard_test.m @@ -268,13 +268,13 @@ function asgard_fokkerplanck1_5p1a_implicit_test(testCase) function asgard_BDF2_test(testCase) addpath(genpath(pwd)); disp('Testing BDF2'); -[err,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... - 'lev',3,'num_steps',2,'dt',0.001,'quiet',true); -[err1,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... - 'lev',3,'num_steps',4,'dt',0.0005,'quiet',true); -[err2,act_f,act_frs] = asgard(diffusion1,'timestep_method','BDF2',... - 'lev',3,'num_steps',8,'dt',0.00025,'quiet',true); -logslope = (log(err1) - log(err))/(log(1/0.001) - log(1/0.0005)); +[err,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... + 'lev',3, 'deg', 4,'num_steps',2,'dt',0.0001,'quiet',true); +[err1,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... + 'lev',3, 'deg', 4,'num_steps',4,'dt',0.00005,'quiet',true); +[err2,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... + 'lev',3, 'deg', 4,'num_steps',8,'dt',0.000025,'quiet',true); +logslope = (log(err1) - log(err))/(log(1/0.0001) - log(1/0.00005)); slopeerr = abs(logslope - 2.000); verifyLessThan(testCase,slopeerr,1e-3); end From 4beafe9c71e2f8308007cebefaf94c9faf854c40 Mon Sep 17 00:00:00 2001 From: lwonnell Date: Thu, 28 May 2020 12:38:51 -0400 Subject: [PATCH 4/4] Added test for BDF2. This test now passes for fokkerplanck4p1_a, showing that it is second-order accurate, but fails for fokkerplanck2_6p1. --- asgard_test.m | 13 +++++-------- time_advance.m | 10 +++++----- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/asgard_test.m b/asgard_test.m index cd982d2e..a9ed5a90 100644 --- a/asgard_test.m +++ b/asgard_test.m @@ -268,14 +268,11 @@ function asgard_fokkerplanck1_5p1a_implicit_test(testCase) function asgard_BDF2_test(testCase) addpath(genpath(pwd)); disp('Testing BDF2'); -[err,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... - 'lev',3, 'deg', 4,'num_steps',2,'dt',0.0001,'quiet',true); -[err1,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... - 'lev',3, 'deg', 4,'num_steps',4,'dt',0.00005,'quiet',true); -[err2,act_f,act_frs] = asgard(advection1,'timestep_method','BDF2',... - 'lev',3, 'deg', 4,'num_steps',8,'dt',0.000025,'quiet',true); -logslope = (log(err1) - log(err))/(log(1/0.0001) - log(1/0.00005)); +[err,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.01,'num_steps',1, 'quiet', true); +[err1,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.005,'num_steps',2,'quiet', true); +[err2,act_f,act_frs] = asgard(fokkerplanck1_4p1a,'lev',4,'deg',5,'timestep_method','BDF2', 'dt',0.0025,'num_steps',4,'quiet', true); +logslope = (log(err1) - log(err))/(log(1/0.1) - log(1/0.05)); slopeerr = abs(logslope - 2.000); -verifyLessThan(testCase,slopeerr,1e-3); +verifyLessThan(testCase,slopeerr,2.5e-2); end diff --git a/time_advance.m b/time_advance.m index f49bba29..01ea9980 100644 --- a/time_advance.m +++ b/time_advance.m @@ -244,11 +244,11 @@ function f2 = bdf2(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax) -s1 = source_vector(pde,opts,hash_table,t+dt); -s2 = source_vector(pde,opts,hash_table,t+2*dt); +s1 = source_vector(pde,opts,hash_table,t+dt/2); +s2 = source_vector(pde,opts,hash_table,t+dt); -bc1 = boundary_condition_vector(pde,opts,hash_table,t+dt); -bc2 = boundary_condition_vector(pde,opts,hash_table,t+2*dt); +bc1 = boundary_condition_vector(pde,opts,hash_table,t+dt/2); +bc2 = boundary_condition_vector(pde,opts,hash_table,t+dt); % %% % Apply any non-identity LHS mass matrix coefficient @@ -256,7 +256,7 @@ applyLHS = ~isempty(pde.termsLHS); %Obtain f1 through same process as Backward Euler -if t <= dt +if t == 0 % opts.timestep_method = 'ode45'; % f1 = ODEm(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax); fac = 200;