From 33c135fc7e225fcaab28cab47cba1090e63e84d4 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Fri, 4 Mar 2022 09:30:40 -0500 Subject: [PATCH 01/12] Reverting IMEX changes from push to reference --- OPTS.m | 2 +- pdes/vlasov_lb_full_f.m | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/OPTS.m b/OPTS.m index 77033d2f..83f08c70 100644 --- a/OPTS.m +++ b/OPTS.m @@ -53,7 +53,7 @@ valid_grid_types = {'SG','FG'}; check_grid_type = @(x) any(validatestring(x,valid_grid_types)); - valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode23s','ode45','RK3','FE','time_independent','matrix_exponential'}; + valid_timestep_methods = {'BE','CN','ode15i','ode15s','ode23s','ode45','RK3','FE','time_independent','matrix_exponential','IMEX'}; check_timestep_method = @(x) any(strcmp(x,valid_timestep_methods)); valid_output_grids = {'quadrature','fixed','uniform','quadrature_with_end_points','dual_valued','elements','interp'}; check_output_grid = @(x) any(strcmp(x,valid_output_grids)); diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index 35d679a4..63a17a6e 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -152,16 +152,16 @@ % div_v(th q) % q = \grad_v f -g1 = @(x,p,t,dat) 1; -g2 = @(x,p,t,dat) 0*x+p.th(x)*p.nu; +g1 = @(x,p,t,dat) 0*x+sqrt(p.th(x)*p.nu); +g2 = @(x,p,t,dat) 0*x+sqrt(p.th(x)*p.nu); pterm1 = MASS(g1,'','',dV); pterm2 = MASS(g2,'','',dV); term_x = SD_TERM({pterm1,pterm2}); g1 = @(x,p,t,dat) 0*x+1; g2 = @(x,p,t,dat) 0*x+1; -pterm1 = DIV(num_dims,g1,'',0,'D','D','','','',dV); -pterm2 = GRAD(num_dims,g2,'',0,'D','D','','','',dV); +pterm1 = DIV(num_dims,g1,'',-1,'D','D','','','',dV); +pterm2 = GRAD(num_dims,g2,'',+1,'D','D','','','',dV); term_v = SD_TERM({pterm1,pterm2}); term5 = MD_TERM(num_dims,{term_x,term_v},'I'); From 8bee45c3a4a7989a81a7ac9670381fe178bb5aef Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Tue, 29 Mar 2022 13:46:54 -0400 Subject: [PATCH 02/12] Got adaptivity working --- fast_2d_matrix_apply.m | 11 ++++-- hash_table_2D_to_1D.m | 55 ++++++++++++++++++------------ moment_reduced_matrix.m | 4 +-- pdes/vlasov_lb_full_f.m | 74 +++++++++++++++++++++++++---------------- pterm_coeff_matrix.m | 6 ++++ sg_to_fg_mapping_2d.m | 6 ++-- time_advance.m | 52 +++++++++++++++++++++-------- 7 files changed, 137 insertions(+), 71 deletions(-) diff --git a/fast_2d_matrix_apply.m b/fast_2d_matrix_apply.m index 5b9c8b41..194731e6 100644 --- a/fast_2d_matrix_apply.m +++ b/fast_2d_matrix_apply.m @@ -9,6 +9,9 @@ if isempty(perm) [perm,iperm,pvec] = sg_to_fg_mapping_2d(pde,opts,A_data); end +if opts.adapt + [perm,iperm,pvec] = sg_to_fg_mapping_2d(pde,opts,A_data); +end % if strcmp(opts.timestep_method,'IMEX') @@ -24,8 +27,10 @@ end %Get dimensions -n = size(pde.terms{1}.terms_1D{1}.mat,1); -m = size(pde.terms{1}.terms_1D{2}.mat,1); +%n = size(pde.terms{1}.terms_1D{1}.mat,1); +%m = size(pde.terms{1}.terms_1D{2}.mat,1); +n = 2^pde.dimensions{1}.lev*opts.deg; +m = 2^pde.dimensions{2}.lev*opts.deg; %Convert to standard FG index f0_F = zeros(numel(perm),1); @@ -37,7 +42,7 @@ %Evaluate sum of krons f1_M = zeros(size(f0_M)); for i=1:numel(term_idx) - f1_M = f1_M + pde.terms{term_idx(i)}.terms_1D{2}.mat*f0_M*pde.terms{term_idx(i)}.terms_1D{1}.mat'; + f1_M = f1_M + pde.terms{term_idx(i)}.terms_1D{2}.mat(1:m,1:m)*f0_M*pde.terms{term_idx(i)}.terms_1D{1}.mat(1:n,1:n)'; end %Convert to vector diff --git a/hash_table_2D_to_1D.m b/hash_table_2D_to_1D.m index 07e73139..439214b5 100644 --- a/hash_table_2D_to_1D.m +++ b/hash_table_2D_to_1D.m @@ -1,9 +1,10 @@ function [hash_table_1D] = hash_table_2D_to_1D(hash_table,opts) -%Converts a 2D hash table into a 1D hash table. This new hash table will -%be a full grid space, but this will eventually be modified to make a 4D -%table into a 2D table. +%Converts a 2D hash table into a 1D full-grid hash_table with standard +%kronecker indexing -assert(size(hash_table.elements.lev_p1,2) == 2); +% TODO -- Actual 2D to 1D hash table conversion + +%assert(size(hash_table.elements.lev_p1,2) == 2); %num_basis = numel(hash_table.elements_idx); %x_dim = 1; @@ -16,26 +17,26 @@ %with respect to that dimension's full grid indexing. % ( This is the same as A_data but % is put here for sanity's sake and reference later ) -%ele_data = full(hash_table.elements.lev_p1(hash_table.elements_idx,:) + ... -% hash_table.elements.pos_p1(hash_table.elements_idx,:) - 1); - -%Create hash_table_1D element_idx based on the maximum element index in the -%x direction of the 2D hash_table. -%hash_table_1D.elements_idx = nan(1,max(ele_data(:,x_dim))); - -%Get maxiumum level for sparse matrix. -max_lev = opts.max_lev; - -%Create sparse containers -%hash_table_1D.elements.lev_p1 = sparse([],[],[],2^max_lev,1,numel(hash_table_1D.elements_idx)); -%hash_table_1D.elements.pos_p1 = sparse([],[],[],2^max_lev,1,numel(hash_table_1D.elements_idx)); - +% ele_data = full(2.^(hash_table.elements.lev_p1(hash_table.elements_idx,:)) + + ... +% hash_table.elements.pos_p1(hash_table.elements_idx,:) - 1); +% +% %Create hash_table_1D element_idx based on the maximum element index in the +% %x direction of the 2D hash_table. +% hash_table_1D.elements_idx = nan(1,max(ele_data(:,x_dim))); +% +% %Get maxiumum level for sparse matrix. +% max_lev = opts.max_lev; +% +% %Create sparse containers +% hash_table_1D.elements.lev_p1 = sparse([],[],[],2^max_lev,1,numel(hash_table_1D.elements_idx)); +% hash_table_1D.elements.pos_p1 = sparse([],[],[],2^max_lev,1,numel(hash_table_1D.elements_idx)); +% % %Populate % count = 1; % for i=1:num_basis % idx = hash_table.elements_idx(i); % if idx <= 2^max_lev -% hash_table_1D.elements_idx(count) = count; %Ordering is standard +% hash_table_1D.elements_idx(count) = idx; % hash_table_1D.elements.lev_p1(count,1) = hash_table.elements.lev_p1(idx,1); % hash_table_1D.elements.pos_p1(count,1) = hash_table.elements.pos_p1(idx,1); % count = count + 1; @@ -43,9 +44,19 @@ % end % hash_table_1D.elements_idx(count:end) = []; -hash_table_1D.elements.lev_p1 = hash_table.elements.lev_p1(1:2^max_lev,1); -hash_table_1D.elements.pos_p1 = hash_table.elements.pos_p1(1:2^max_lev,1); -hash_table_1D.elements_idx = 1:nnz(hash_table_1D.elements.pos_p1); +%Get maximum level for x +max_lev = max(hash_table.elements.lev_p1(:,1))-1; + +%Create full grid standard hash_table with max_lev +hash_table_1D.elements_idx = 1:2^(max_lev); +hash_table_1D.elements.lev_p1 = sparse(ceil(log2(1:2^max_lev)))'+1; +hash_table_1D.elements.pos_p1 = (1:2^max_lev)' - ... + floor(2.^(hash_table_1D.elements.lev_p1-2)); %Time for some magic... + + +% hash_table_1D.elements.lev_p1 = hash_table.elements.lev_p1(1:2^max_lev,1); +% hash_table_1D.elements.pos_p1 = hash_table.elements.pos_p1(1:2^max_lev,1); +% hash_table_1D.elements_idx = 1:nnz(hash_table_1D.elements.pos_p1); end diff --git a/moment_reduced_matrix.m b/moment_reduced_matrix.m index d6c3e76c..da67b8a3 100644 --- a/moment_reduced_matrix.m +++ b/moment_reduced_matrix.m @@ -23,8 +23,8 @@ -%Get row-dim of M -rows = max(A_data.element_local_index_D{x_dim}*deg); +%Get row-dim of M -- Assuming a full grid output here +rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); %Containers for sparse matrix I = zeros(deg^2*num_ele,1); %Square is hardcoded. Needs to change diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index 63a17a6e..6b324e55 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -13,10 +13,6 @@ % asgard(@vlasov_lb_full_f,'timestep_method','IMEX','deg',3,'lev',[8 3],'dt',0.0002,'num_steps',500,'grid_type','FG','output_grid','interp','quiet',true,'build_realspace_output',false) % -soln_x = @(x,p,t) 0*x+1; -soln_v = @(v,p,t) 0*v+1; -soln_t = @(t,p) 0*t+1; - %% Define the dimensions % % Here we setup a 2D problem (x,v) @@ -52,14 +48,31 @@ %% Construct (n,u,theta) -params.n = @(x) 1*(x<-0.5) + 1/8*(x >= -0.5).*(x <= 0.5) + 1*(x>0.5); -params.u = @(x) 0; -params.th = @(x) 1*(x<-0.5) + 4/5*(x >= -0.5).*(x <= 0.5) + 1*(x>0.5); - -%params.n = @(x) 0*x+1; -%params.u = @(x) 0*x; -%params.th = @(x) 0*x+1; -params.nu = 1e3; +switch opts.case_ + case 1 + params.n = @(x) 1*(x<-0.5) + 1/8*(x >= -0.5).*(x <= 0.5) + 1*(x>0.5); + params.u = @(x) 0; + params.th = @(x) 1*(x<-0.5) + 4/5*(x >= -0.5).*(x <= 0.5) + 1*(x>0.5); + params.nu = 1e3; + case 2 + % deg = 3, lev = [8,5] + U_0 = 1; + T_0 = 1e-0; + + params.n = @(x) 0.5*sin(2*pi*x)+1; + params.u = @(x) 0*x+U_0; + params.th = @(x) 0*x+T_0; + params.nu = 0; + case 3 + N_0 = 1; + T_0 = 1/3; + C_0 = sqrt(3*T_0); + Amp = 1e-6; + + params.n = @(x) N_0 + Amp*sin(2*pi*x)/C_0^2; + params.u = @(x) Amp*sin(2*pi*x)./(C_0*params.n(x)); + params.th = @(x) (N_0*T_0+Amp*sin(2*pi*x))./params.n(x)-params.u(x).^2; +end %% Define the analvtic solution (optional). @@ -69,18 +82,23 @@ solutions = {}; %% Initial conditions -ic1 = new_md_func(num_dims,{@(x,p,t) (abs(x) > 0.5),... - @(v,p,t) 1/sqrt(2*pi)*exp(-v.^2/2),... - @(t,p) 0*t+1}); -ic2 = new_md_func(num_dims,{@(x,p,t) (abs(x) <= 0.5),... - @(v,p,t) (1/8)/sqrt(2*pi*4/5)*exp(-v.^2/(2*4/5)),... - @(t,p) 0*t+1}); -initial_conditions = {ic1,ic2}; - -%ic1 = new_md_func(num_dims,{@(x,p,t) 0*x+1,... -% @(v,p,t) 1/sqrt(2*pi)*exp(-v.^2/2),... -% @(t,p) 0*t+1}); -%initial_conditions = {ic1}; + +switch opts.case_ + case 1 + ic1 = new_md_func(num_dims,{@(x,p,t) (abs(x) > 0.5),... + @(v,p,t) 1/sqrt(2*pi)*exp(-v.^2/2),... + @(t,p) 0*t+1}); + ic2 = new_md_func(num_dims,{@(x,p,t) (abs(x) <= 0.5),... + @(v,p,t) (1/8)/sqrt(2*pi*4/5)*exp(-v.^2/(2*4/5)),... + @(t,p) 0*t+1}); + initial_conditions = {ic1,ic2}; + case 2 + %U_0 and T_0 defined above + ic1 = new_md_func(num_dims,{@(x,p,t) p.n(x)/sqrt(2*pi*T_0),... + @(v,p,t) exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + initial_conditions = {ic1}; +end %% Define the terms of the PDE % @@ -152,16 +170,16 @@ % div_v(th q) % q = \grad_v f -g1 = @(x,p,t,dat) 0*x+sqrt(p.th(x)*p.nu); -g2 = @(x,p,t,dat) 0*x+sqrt(p.th(x)*p.nu); +g1 = @(x,p,t,dat) 0*x+1; +g2 = @(x,p,t,dat) 0*x+p.th(x)*p.nu; pterm1 = MASS(g1,'','',dV); pterm2 = MASS(g2,'','',dV); term_x = SD_TERM({pterm1,pterm2}); g1 = @(x,p,t,dat) 0*x+1; g2 = @(x,p,t,dat) 0*x+1; -pterm1 = DIV(num_dims,g1,'',-1,'D','D','','','',dV); -pterm2 = GRAD(num_dims,g2,'',+1,'D','D','','','',dV); +pterm1 = DIV(num_dims,g1,'',0,'D','D','','','',dV); +pterm2 = GRAD(num_dims,g2,'',0,'D','D','','','',dV); term_v = SD_TERM({pterm1,pterm2}); term5 = MD_TERM(num_dims,{term_x,term_v},'I'); diff --git a/pterm_coeff_matrix.m b/pterm_coeff_matrix.m index 5b50316f..d60fd2a8 100644 --- a/pterm_coeff_matrix.m +++ b/pterm_coeff_matrix.m @@ -13,6 +13,12 @@ [M,M_not_rotated] = coeff_matrix(deg,t,dim,lhs_mass_pterm,params,FMWT_blocks,coeff_level); pterm.LHS_mass_mat = full(M); pterm.LHS_mass_mat_unrotated = full(M_not_rotated); + identity_ref = norm(pterm.LHS_mass_mat-eye(size(M)),'fro'); + if identity_ref < 1e-11 %%Checking for identity mass matrix + %fprintf('Replacing mass mat with identity. Error = %e\n',identity_ref); + pterm.LHS_mass_mat = speye(size(M)); + pterm.LHS_mass_mat_unrotated = speye(size(M)); + end end % Move M to the RHS diff --git a/sg_to_fg_mapping_2d.m b/sg_to_fg_mapping_2d.m index 62602cbd..2b9578c2 100644 --- a/sg_to_fg_mapping_2d.m +++ b/sg_to_fg_mapping_2d.m @@ -46,8 +46,10 @@ iperm = zeros(num_sparse_ele*deg^2,1); %Get FG_dof for each direction -max_x_dof = size(pde.terms{1}.terms_1D{1}.mat,1); -max_y_dof = size(pde.terms{1}.terms_1D{2}.mat,1); +%max_x_dof = size(pde.terms{1}.terms_1D{1}.mat,1); +%max_y_dof = size(pde.terms{1}.terms_1D{2}.mat,1); +max_x_dof = 2^pde.dimensions{1}.lev*deg; +max_y_dof = 2^pde.dimensions{2}.lev*deg; perm = nan(max_x_dof*max_y_dof,1); idx = 1; %Keep track of index diff --git a/time_advance.m b/time_advance.m index 41b7cc28..4668a17d 100644 --- a/time_advance.m +++ b/time_advance.m @@ -498,10 +498,6 @@ function f1 = imex(pde,opts,A_data,f0,t,dt,deg,hash_table,Vmax,Emax) -persistent pde_1d -persistent nodes -persistent Meval - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% C++ Implementation Notes %%%%%%%%%%%%%%%%%%%%%%% @@ -516,6 +512,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +persistent P %Switch for IMEX iteration: % 1 : Backward Euler in explicit terms 'E', @@ -526,16 +523,12 @@ BEFE = 0; -if isempty(Meval) - %Compute everything that will not change per time iteration +%Get quadrature points in realspace stiffness matrix calculation +[Meval,nodes] = matrix_plot_D(pde,opts,pde.dimensions{1}); + +%Make fake pde file for use in physical realspace transformation +pde_1d.dimensions = pde.dimensions(1); - %Get quadrature points in realspace stiffness matrix calculation - [Meval,nodes] = matrix_plot_D(pde,opts,pde.dimensions{1}); - - %Make fake pde file for use in physical realspace transformation - pde_1d.dimensions = pde.dimensions(1); - -end %Create moment matrices that take DG function in (x,v) and transfer it %to DG function in x. @@ -546,6 +539,12 @@ end end +%if isempty(P) +% M = [moment_mat{1};moment_mat{2};moment_mat{3}]; +% [~,S,V] = svd(full(M),'econ'); +% P = V(:,1:sum(diag(S) > 1e-12)); +%end + hash_table_1D = hash_table_2D_to_1D(hash_table,opts); assert(isempty(pde.termsLHS),'LHS terms currently not supported by IMEX'); @@ -612,7 +611,8 @@ %%%%% %Explicit step - f_2s = f0 + dt*fast_2d_matrix_apply(opts,pde,A_data,f0,'E'); + f_2s = f0 + dt*fast_2d_matrix_apply(opts,pde,A_data,f0,'E'); + %f_2s = f0; %Create rho_2s mom0 = moment_mat{1}*f_2s; %integral of (f,1)_v @@ -626,6 +626,26 @@ mom2 = moment_mat{3}*f_2s; %integral of (f,v^2)_v mom2_real = wavelet_to_realspace(pde_1d,opts,{Meval},mom2,hash_table_1D); pde.params.th = @(x) interp1(nodes,mom2_real,x,'nearest','extrap')./pde.params.n(x) - pde.params.u(x).^2; + + %Plot moments + if ~opts.quiet || 1 + + fig1 = figure(1000); + fig1.Units = 'Normalized'; + fig1.Position = [0.5 0.5 0.3 0.3]; + subplot(2,2,1); + plot(nodes,pde.params.n(nodes)); + title('n_f'); + subplot(2,2,2); + plot(nodes,pde.params.u(nodes)); + title('u_f'); + subplot(2,2,3); + plot(nodes,pde.params.th(nodes)); + title('th_f'); + sgtitle("Fluid Variables. t = "+num2str(t+dt)); + drawnow + + end %Update coefficients pde = get_coeff_mats(pde,opts,t,0); @@ -636,6 +656,7 @@ fprintf('BICGSTABL did not converge. flag = %d, relres = %5.4e\n',flag,relres); assert(relres < 1e-10) end + %f_2 = f_2 - P*(P'*(f_2-f_2s)); %%%%% %%% Third stage @@ -666,6 +687,7 @@ fprintf('BICGSTABL did not converge. flag = %d, relres = %5.4e\n',flag,relres); assert(relres < 1e-10) end + %f_3 = f_3 - P*(P'*(f_3-f_3s)); %Update timestep to final stage f1 = f_3; @@ -690,6 +712,8 @@ end + %fprintf("Norm Check t = %f, |f1|_2 = %f\n",t+dt,norm(f1)); + end end From 0c1f36a864eaf5e451297c75951eadd8e6695400 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Tue, 29 Mar 2022 17:06:04 -0400 Subject: [PATCH 03/12] Made moment_reduced_matrix work for 1x-2v --- moment_reduced_matrix.m | 106 ++++++++++++++++++++++++++++++---------- 1 file changed, 79 insertions(+), 27 deletions(-) diff --git a/moment_reduced_matrix.m b/moment_reduced_matrix.m index da67b8a3..9c20a317 100644 --- a/moment_reduced_matrix.m +++ b/moment_reduced_matrix.m @@ -1,9 +1,9 @@ function [M] = moment_reduced_matrix(opts,pde,A_data,hash_table,moment_idx) % % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% CURRENTLY ONLY WORKS FOR 1x-1v system -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% CURRENTLY ONLY WORKS FOR 1x-1v or 1x-2v system +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % This function produces a matrix M that is the transformation matrix that @@ -15,41 +15,93 @@ % In 1x-1v systems, the output function is full grid, but this will break % down in multiple dimensions. +num_dims = numel(pde.dimensions); deg = opts.deg; num_ele = numel(A_data.element_local_index_D{1}); -x_dim = 1; %Hardcoded for now. Needs to change in multi-dim -v_dim = 2; %Hardcoded for now. Needs to change in multi-dim -g_vec = pde.moments{moment_idx}.fList{v_dim}; +if num_dims == 2 % 1x-1v + x_dim = 1; %Hardcoded for now. Needs to change in multi-dim + v_dim = 2; %Hardcoded for now. Needs to change in multi-dim + g_vec = pde.moments{moment_idx}.fList{v_dim}; -%Get row-dim of M -- Assuming a full grid output here -rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); + %Get row-dim of M -- Assuming a full grid output here + rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); -%Containers for sparse matrix -I = zeros(deg^2*num_ele,1); %Square is hardcoded. Needs to change -J = zeros(deg^2*num_ele,1); -S = zeros(deg^2*num_ele,1); + %Containers for sparse matrix + I = zeros(deg^2*num_ele,1); %Square is hardcoded. Needs to change + J = zeros(deg^2*num_ele,1); + S = zeros(deg^2*num_ele,1); -count = 1; -k = 1:deg; -for i=1:num_ele - for j=1:deg - temp = ((A_data.element_local_index_D{x_dim}(i)-1)*deg + j) + 0*k; - I(count:count+deg-1) = temp'; - - temp = (i-1)*deg^2 + (j-1)*deg + k; - J(count:count+deg-1) = temp'; - - temp = g_vec((A_data.element_local_index_D{v_dim}(i)-1)*deg + k); - S(count:count+deg-1) = temp'; + count = 1; + k = 1:deg; %Not going to loop over degree + for i=1:num_ele + for j=1:deg + temp = ((A_data.element_local_index_D{x_dim}(i)-1)*deg + j) + 0*k; + I(count:count+deg-1) = temp'; + + temp = (i-1)*deg^2 + (j-1)*deg + k; + J(count:count+deg-1) = temp'; + + temp = g_vec((A_data.element_local_index_D{v_dim}(i)-1)*deg + k); + S(count:count+deg-1) = temp'; + + count = count + deg; + end + end + + M = sparse(I,J,S,rows,deg^2*num_ele); + +elseif num_dims == 3 % 1x-2v + + x_dim = 1; %Hardcoded for now. Needs to change in multi-dim + v_dim_1 = 2; %Hardcoded for now. Needs to change in multi-dim + v_dim_2 = 2; %Hardcoded for now. Needs to change in multi-dim + + %Get row-dim of M -- Assuming a full grid output here + rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); + + %Containers for sparse matrix + I = zeros(deg^3*num_ele,1); %Square is hardcoded. Needs to change + J = zeros(deg^3*num_ele,1); + S = zeros(deg^3*num_ele,1); + + g_vec_1 = pde.moments{moment_idx}.fList{v_dim_1}; + g_vec_2 = pde.moments{moment_idx}.fList{v_dim_2}; + + count = 0; + for i=1:num_ele %loop through all elements lit up - count = count + deg; + %Get local dof for element + l_dof = [A_data.element_local_index_D{x_dim }(i)-1,... + A_data.element_local_index_D{v_dim_1}(i)-1,... + A_data.element_local_index_D{v_dim_2}(i)-1]; + + for xdeg=1:deg + %The row of the matrix is given by the x-dof + row_idx = l_dof(1)*deg+xdeg; + for vdeg1=1:deg + for vdeg2=1:deg + %Get column entry based on degree index + col_idx = (i-1)*deg^3 + (xdeg-1)*deg^2 + (vdeg1-1)*deg + vdeg2; + %Get entry (which is the product of the 1d elements of + %g_vec) + val = g_vec_1((A_data.element_local_index_D{v_dim_1}(i)-1)*deg+vdeg1) * ... + g_vec_2((A_data.element_local_index_D{v_dim_2}(i)-1)*deg+vdeg2); + I(count+1) = row_idx; + J(count+1) = col_idx; + S(count+1) = val; + count = count + 1; + end + end + end end + + %Construct sparse matrix + M = sparse(I,J,S,rows,deg^3*num_ele); + end -M = sparse(I,J,S,rows,deg^2*num_ele); - end From d7c3def30a9c47b8b7e9dc7b5fe7c2af5283c6d9 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Thu, 31 Mar 2022 08:59:04 -0400 Subject: [PATCH 04/12] Added analytic soln to case 2 of vlasov_full_f --- pdes/vlasov_lb_full_f.m | 19 ++++++++++- time_advance.m | 76 ++++++++++++++++++++++++++++++----------- 2 files changed, 74 insertions(+), 21 deletions(-) diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index 6b324e55..cbbee81d 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -79,7 +79,24 @@ %soln1 = new_md_func(num_dims,{soln_x,soln_v,soln_t}); %solutions = {soln1}; -solutions = {}; +switch opts.case_ + case 2 + soln1 = new_md_func(num_dims,... + {@(x,p,t,dat) 0.5*sin(2*pi*x),... + @(v,p,t,dat) cos(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + soln2 = new_md_func(num_dims,... + {@(x,p,t,dat) -0.5*cos(2*pi*x),... + @(v,p,t,dat) sin(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + soln3 = new_md_func(num_dims,... + {@(x,p,t,dat) 0*x+1,... + @(v,p,t,dat) 1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + solutions = {soln1,soln2,soln3}; + otherwise + solutions = {}; +end %% Initial conditions diff --git a/time_advance.m b/time_advance.m index 4668a17d..cf4e214e 100644 --- a/time_advance.m +++ b/time_advance.m @@ -628,24 +628,24 @@ pde.params.th = @(x) interp1(nodes,mom2_real,x,'nearest','extrap')./pde.params.n(x) - pde.params.u(x).^2; %Plot moments - if ~opts.quiet || 1 - - fig1 = figure(1000); - fig1.Units = 'Normalized'; - fig1.Position = [0.5 0.5 0.3 0.3]; - subplot(2,2,1); - plot(nodes,pde.params.n(nodes)); - title('n_f'); - subplot(2,2,2); - plot(nodes,pde.params.u(nodes)); - title('u_f'); - subplot(2,2,3); - plot(nodes,pde.params.th(nodes)); - title('th_f'); - sgtitle("Fluid Variables. t = "+num2str(t+dt)); - drawnow - - end +% if ~opts.quiet || 1 +% +% fig1 = figure(1000); +% fig1.Units = 'Normalized'; +% fig1.Position = [0.5 0.5 0.3 0.3]; +% subplot(2,2,1); +% plot(nodes,pde.params.n(nodes)); +% title('n_f'); +% subplot(2,2,2); +% plot(nodes,pde.params.u(nodes)); +% title('u_f'); +% subplot(2,2,3); +% plot(nodes,pde.params.th(nodes)); +% title('th_f'); +% sgtitle("Fluid Variables. t = "+num2str(t+dt)); +% drawnow +% +% end %Update coefficients pde = get_coeff_mats(pde,opts,t,0); @@ -692,8 +692,44 @@ %Update timestep to final stage f1 = f_3; - %Plot moments - if ~opts.quiet || 1 + + if opts.case_ == 2 + %Get analytic solution at t_{n+1} + fval_analytic = exact_solution_vector(pde,opts,hash_table,t+dt); + + %Get analytic moments + + %n + mom0_a = moment_mat{1}*fval_analytic; + mom0_a_r = wavelet_to_realspace(pde_1d,opts,{Meval},mom0_a,hash_table_1D); + analytic_moments.n = mom0_a_r; + + %u + mom1_a = moment_mat{2}*fval_analytic; + mom1_a_r = wavelet_to_realspace(pde_1d,opts,{Meval},mom1_a,hash_table_1D); + analytic_moments.u = mom1_a_r./mom0_a_r; + + %T + mom2_a = moment_mat{3}*fval_analytic; + mom2_a_r = wavelet_to_realspace(pde_1d,opts,{Meval},mom2_a,hash_table_1D); + analytic_moments.T = mom2_a_r./mom0_a_r - (analytic_moments.u).^2; + + fig1 = figure(1000); + fig1.Units = 'Normalized'; + fig1.Position = [0.5 0.5 0.3 0.3]; + subplot(2,2,1); + plot(nodes,pde.params.n(nodes),nodes,analytic_moments.n); + title('n_f'); + subplot(2,2,2); + plot(nodes,pde.params.u(nodes),nodes,analytic_moments.u); + title('u_f'); + subplot(2,2,3); + plot(nodes,pde.params.th(nodes),nodes,analytic_moments.T); + title('th_f'); + sgtitle("Fluid Variables. t = "+num2str(t+dt)); + drawnow + else + %Plot moments fig1 = figure(1000); fig1.Units = 'Normalized'; From b386a7045e08aa6e6f5ea091117b508de2e854b5 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Tue, 5 Apr 2022 16:27:29 -0400 Subject: [PATCH 05/12] Adding moment_test for C++ conversion --- moment_test.m | 100 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 moment_test.m diff --git a/moment_test.m b/moment_test.m new file mode 100644 index 00000000..c9462e32 --- /dev/null +++ b/moment_test.m @@ -0,0 +1,100 @@ +function tests = moment_test() +tests = functiontests(localfunctions); +end + +function test_moment_2D_spherical(testCase) + +args = {'lev',5,'deg',4,'case',3}; +opts = OPTS(args); +pde = diffusion2(opts); + +num_dims = numel(pde.dimensions); +[elements, elements_idx] = hash_table_sparse_nD (pde.get_lev_vec, opts.max_lev, opts.grid_type); +hash_table.elements = elements; +hash_table.elements_idx = elements_idx; +t = 0; + +%% Verify moment f(x,y) = 1 (should be auto added) + +ic_x = @(x,p,t) cos(pi*x); +ic_y = @(y,p,t) cos(pi*y); +ic_t = @(t,p) exp(-2*pi^2*t); +pde.initial_conditions = {new_md_func(num_dims,{ic_x,ic_y,ic_t})}; + +pde = compute_dimension_mass_mat(opts,pde); +pde = calculate_moment_data(pde,opts); + +fval = initial_condition_vector(pde,opts,hash_table,t); +[pde, moment_val, ~] = calculate_mass(pde,opts,fval,hash_table,t); + +%moment_val should be 0 +diff1 = abs(moment_val-0); +verifyLessThan(testCase, diff1, 1e-10); + +%% Verify moment f(x,y) = ic + +pde.moments{1} = MOMENT({new_md_func(num_dims,{ic_x,ic_y,ic_t})}); + +pde = compute_dimension_mass_mat(opts,pde); +pde = calculate_moment_data(pde,opts); + +fval = initial_condition_vector(pde,opts,hash_table,t); +[pde, moment_val, ~] = calculate_mass(pde,opts,fval,hash_table,t); + +%moment_val should be 0 +diff1 = abs(moment_val-0.25); +verifyLessThan(testCase, diff1, 1e-13); + +%% Verify moment f(x,y) = y + +ic_x = @(x,p,t) 1-x.^2; +ic_y = @(y,p,t) y; +ic_t = @(t,p) 0*t+1; +pde.initial_conditions = {new_md_func(num_dims,{ic_x,ic_y,ic_t})}; + +mom_x = @(x,p,t) 0*x+1; +mom_y = @(y,p,t) y; +mom_t = @(t,p) 0*t+1; +mom_func = new_md_func(num_dims,{mom_x,mom_y,mom_t}); + +pde.moments{1} = MOMENT({mom_func}); + +pde = compute_dimension_mass_mat(opts,pde); +pde = calculate_moment_data(pde,opts); + +fval = initial_condition_vector(pde,opts,hash_table,t); +[pde, moment_val, ~] = calculate_mass(pde,opts,fval,hash_table,t); + +%moment_val should be 2/9 +diff1 = abs(moment_val-2/9); +verifyLessThan(testCase, diff1, 1e-13); + + +%% Do non-cartesian coordinate example + +pde = diffusion2_spherical(opts); + +num_dims = numel(pde.dimensions); +[elements, elements_idx] = hash_table_sparse_nD (pde.get_lev_vec, opts.max_lev, opts.grid_type); +hash_table.elements = elements; +hash_table.elements_idx = elements_idx; +t = 0; + +mom_x = @(r,p,t) 0*r+1; +mom_y = @(th,p,t) cos(th); +mom_t = @(t,p) 0*t+1; +mom_func = new_md_func(num_dims,{mom_x,mom_y,mom_t}); +pde.moments{1} = MOMENT({mom_func}); + +pde = compute_dimension_mass_mat(opts,pde); +pde = calculate_moment_data(pde,opts); + +fval = initial_condition_vector(pde,opts,hash_table,t); +[pde, moment_val, ~] = calculate_mass(pde,opts,fval,hash_table,t); + +%moment_val should be 0 +diff1 = abs(moment_val-8/3); +verifyLessThan(testCase, diff1, 1e-13); + + +end \ No newline at end of file From aa2e34b2f1539bca14c4d2162726f70ee7e733ed Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Wed, 4 May 2022 10:27:01 -0400 Subject: [PATCH 06/12] Changes for sparse grid investigation --- pdes/vlasov_lb_full_f.m | 2 +- time_advance.m | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index cbbee81d..f134d276 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -111,7 +111,7 @@ initial_conditions = {ic1,ic2}; case 2 %U_0 and T_0 defined above - ic1 = new_md_func(num_dims,{@(x,p,t) p.n(x)/sqrt(2*pi*T_0),... + ic1 = new_md_func(num_dims,{@(x,p,t) (0.5*sin(2*pi*x)+1)./sqrt(2*pi*T_0),... @(v,p,t) exp(-(v-U_0).^2/(2*T_0)),... @(t,p) 0*t+1}); initial_conditions = {ic1}; diff --git a/time_advance.m b/time_advance.m index cf4e214e..54517490 100644 --- a/time_advance.m +++ b/time_advance.m @@ -664,6 +664,7 @@ f_3s = f0 + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f0+f_2,'E') ... + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f_2,'I'); + %f_3s = f0; %Create rho_3s mom0 = moment_mat{1}*f_3s; %integral of (f,1)_v @@ -693,9 +694,10 @@ f1 = f_3; - if opts.case_ == 2 + if opts.case_ == 2 && 0 %Get analytic solution at t_{n+1} fval_analytic = exact_solution_vector(pde,opts,hash_table,t+dt); + %fval_analytic = exact_solution_vector(pde,opts,hash_table,t); %Get analytic moments @@ -718,13 +720,13 @@ fig1.Units = 'Normalized'; fig1.Position = [0.5 0.5 0.3 0.3]; subplot(2,2,1); - plot(nodes,pde.params.n(nodes),nodes,analytic_moments.n); + plot(nodes,pde.params.n(nodes)-analytic_moments.n'); title('n_f'); subplot(2,2,2); - plot(nodes,pde.params.u(nodes),nodes,analytic_moments.u); + plot(nodes,pde.params.u(nodes)-analytic_moments.u'); title('u_f'); subplot(2,2,3); - plot(nodes,pde.params.th(nodes),nodes,analytic_moments.T); + plot(nodes,pde.params.th(nodes)-analytic_moments.T'); title('th_f'); sgtitle("Fluid Variables. t = "+num2str(t+dt)); drawnow From 6340db81e9b4c6ee1908aba524d1f377c5db6e20 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Thu, 5 May 2022 08:55:10 -0400 Subject: [PATCH 07/12] Only apply realspace transform when plotting --- asgard_run_pde.m | 117 +++++++++++++++++++++++------------------------ time_advance.m | 2 +- 2 files changed, 57 insertions(+), 62 deletions(-) diff --git a/asgard_run_pde.m b/asgard_run_pde.m index 5f9f5de8..22576656 100644 --- a/asgard_run_pde.m +++ b/asgard_run_pde.m @@ -89,36 +89,33 @@ %% Construct transforms back to realspace for plotting -if opts.build_realspace_output - for d=1:num_dims - if strcmp(opts.output_grid,'fixed') - if d==1 - num_fixed_grid = 51; - else - num_fixed_grid = 21; - end - nodes_nodups{d} = ... - linspace(pde.dimensions{d}.min,pde.dimensions{d}.max,num_fixed_grid); - [Meval{d},nodes{d},nodes_count{d}] = ... - matrix_plot_D(pde,opts,pde.dimensions{d},nodes_nodups{d}); - elseif strcmp(opts.output_grid,'elements') - [element_coordinates,element_coordinates_deg] = get_sparse_grid_coordinates(pde,opts,hash_table); - nodes_nodups{d} = unique(sort(element_coordinates_deg(d,:))); - [Meval{d},nodes{d},nodes_count{d}] = ... - matrix_plot_D(pde,opts,pde.dimensions{d},nodes_nodups{d}); +for d=1:num_dims + if strcmp(opts.output_grid,'fixed') + if d==1 + num_fixed_grid = 51; else - [Meval{d},nodes{d}] = matrix_plot_D(pde,opts,pde.dimensions{d}); - nodes_nodups{d} = nodes{d}; - nodes_count{d} = nodes{d}.*0+1; + num_fixed_grid = 21; end + nodes_nodups{d} = ... + linspace(pde.dimensions{d}.min,pde.dimensions{d}.max,num_fixed_grid); + [Meval{d},nodes{d},nodes_count{d}] = ... + matrix_plot_D(pde,opts,pde.dimensions{d},nodes_nodups{d}); + elseif strcmp(opts.output_grid,'elements') + [element_coordinates,element_coordinates_deg] = get_sparse_grid_coordinates(pde,opts,hash_table); + nodes_nodups{d} = unique(sort(element_coordinates_deg(d,:))); + [Meval{d},nodes{d},nodes_count{d}] = ... + matrix_plot_D(pde,opts,pde.dimensions{d},nodes_nodups{d}); + else + [Meval{d},nodes{d}] = matrix_plot_D(pde,opts,pde.dimensions{d}); + nodes_nodups{d} = nodes{d}; + nodes_count{d} = nodes{d}.*0+1; end end + %% Construct a n-D coordinate array -if opts.build_realspace_output - coord = get_realspace_coords(pde,nodes); - coord_nodups = get_realspace_coords(pde,nodes_nodups); -end +coord = get_realspace_coords(pde,nodes); +coord_nodups = get_realspace_coords(pde,nodes_nodups); %% Plot initial condition if num_dims <=3 @@ -126,13 +123,12 @@ %% % Get the real space solution - if opts.build_realspace_output - fval_realspace = wavelet_to_realspace(pde,opts,Meval,fval,hash_table); - if ~isempty(pde.solutions) - fval_realspace_analytic = get_analytic_realspace_solution_D(pde,opts,coord,t); - fval_realspace_analytic = reshape(fval_realspace_analytic, length(fval_realspace),1); - end + fval_realspace = wavelet_to_realspace(pde,opts,Meval,fval,hash_table); + if ~isempty(pde.solutions) + fval_realspace_analytic = get_analytic_realspace_solution_D(pde,opts,coord,t); + fval_realspace_analytic = reshape(fval_realspace_analytic, length(fval_realspace),1); end + % construct the moment function handle list for calculating the mass if opts.calculate_mass @@ -146,40 +142,38 @@ fval_realspace_analytic = reshape(fval_realspace_analytic, length(fval_realspace), 1); end - if opts.build_realspace_output - f_realspace_nD = singleD_to_multiD(num_dims,fval_realspace,nodes); - if strcmp(opts.output_grid,'fixed') || strcmp(opts.output_grid,'elements') - f_realspace_nD = ... - remove_duplicates(num_dims,f_realspace_nD,nodes_nodups,nodes_count); - end + f_realspace_nD = singleD_to_multiD(num_dims,fval_realspace,nodes); + if strcmp(opts.output_grid,'fixed') || strcmp(opts.output_grid,'elements') + f_realspace_nD = ... + remove_duplicates(num_dims,f_realspace_nD,nodes_nodups,nodes_count); + end - if isempty(pde.solutions) - f_realspace_analytic_nD = []; - else - f_realspace_analytic_nD = get_analytic_realspace_solution_D(pde,opts,coord_nodups,t); - end - - if opts.save_output - if num_dims <= 3 - f_realspace_nD_t{1} = f_realspace_nD; - else - error('Save output for num_dimensions >3 not yet implemented'); - end - end + if isempty(pde.solutions) + f_realspace_analytic_nD = []; + else + f_realspace_analytic_nD = get_analytic_realspace_solution_D(pde,opts,coord_nodups,t); + end - if opts.use_oldhash + if opts.save_output + if num_dims <= 3 + f_realspace_nD_t{1} = f_realspace_nD; else - element_coordinates = get_sparse_grid_coordinates(pde,opts,hash_table); + error('Save output for num_dimensions >3 not yet implemented'); end + end - if norm(fval_realspace) > 0 && ~opts.quiet - figs.ic = figure('Name','Initial Condition','Units','normalized','Position',[0.7,0.1,0.3,0.3]); + if opts.use_oldhash + else + element_coordinates = get_sparse_grid_coordinates(pde,opts,hash_table); + end - if opts.use_oldhash - plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD); - else - plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD,element_coordinates); - end + if ~opts.quiet + figs.ic = figure('Name','Initial Condition','Units','normalized','Position',[0.7,0.1,0.3,0.3]); + + if opts.use_oldhash + plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD); + else + plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD,element_coordinates); end end @@ -407,7 +401,7 @@ % Write the present fval to file. if write_fval; write_fval_to_file(fval,lev,deg,L); end - if num_dims <=3 && opts.build_realspace_output + if num_dims <=3 && ( opts.build_realspace_output || (mod(L,opts.plot_freq) == 0) ) %% % Get the real space solution @@ -492,7 +486,7 @@ % Reshape realspace solution and plot - if num_dims <= 3 && opts.build_realspace_output + if num_dims <= 3 && (opts.build_realspace_output || (mod(L,opts.plot_freq) == 0) ) f_realspace_nD = singleD_to_multiD(num_dims,fval_realspace,nodes); if strcmp(opts.output_grid,'fixed') || strcmp(opts.output_grid,'elements') @@ -520,7 +514,8 @@ figs.solution = figure('Name','Solution','Units','normalized','Position',[0.1,0.1,0.3,0.3]); end - plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD,element_coordinates); + %plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD,element_coordinates); + plot_fval(pde,nodes_nodups,f_realspace_nD,f_realspace_analytic_nD); % this is just for the RE paper plot_fval_in_cyl = false; diff --git a/time_advance.m b/time_advance.m index 54517490..2f0a6e2b 100644 --- a/time_advance.m +++ b/time_advance.m @@ -694,7 +694,7 @@ f1 = f_3; - if opts.case_ == 2 && 0 + if opts.case_ == 2 %Get analytic solution at t_{n+1} fval_analytic = exact_solution_vector(pde,opts,hash_table,t+dt); %fval_analytic = exact_solution_vector(pde,opts,hash_table,t); From b6f381ac158fe48c60a39391feda6758b82fb2c7 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Fri, 6 May 2022 12:27:08 -0400 Subject: [PATCH 08/12] changing data for plots --- pdes/vlasov_lb_full_f.m | 32 ++++++++++++++++++-------------- time_advance.m | 11 ++++++----- 2 files changed, 24 insertions(+), 19 deletions(-) diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index f134d276..741c5e28 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -80,20 +80,24 @@ %solutions = {soln1}; switch opts.case_ - case 2 - soln1 = new_md_func(num_dims,... - {@(x,p,t,dat) 0.5*sin(2*pi*x),... - @(v,p,t,dat) cos(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... - @(t,p) 0*t+1}); - soln2 = new_md_func(num_dims,... - {@(x,p,t,dat) -0.5*cos(2*pi*x),... - @(v,p,t,dat) sin(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... - @(t,p) 0*t+1}); - soln3 = new_md_func(num_dims,... - {@(x,p,t,dat) 0*x+1,... - @(v,p,t,dat) 1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... - @(t,p) 0*t+1}); - solutions = {soln1,soln2,soln3}; + case 2 + if params.nu < 1e-8 + soln1 = new_md_func(num_dims,... + {@(x,p,t,dat) 0.5*sin(2*pi*x),... + @(v,p,t,dat) cos(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + soln2 = new_md_func(num_dims,... + {@(x,p,t,dat) -0.5*cos(2*pi*x),... + @(v,p,t,dat) sin(2*pi*v*t).*1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + soln3 = new_md_func(num_dims,... + {@(x,p,t,dat) 0*x+1,... + @(v,p,t,dat) 1/sqrt(2*pi*T_0).*exp(-(v-U_0).^2/(2*T_0)),... + @(t,p) 0*t+1}); + solutions = {soln1,soln2,soln3}; + else + solutions = {}; + end otherwise solutions = {}; end diff --git a/time_advance.m b/time_advance.m index 2f0a6e2b..8f4238e4 100644 --- a/time_advance.m +++ b/time_advance.m @@ -693,10 +693,11 @@ %Update timestep to final stage f1 = f_3; - - if opts.case_ == 2 + time = t+dt; + if opts.case_ == 2 && (abs(pde.params.nu) < 1e-8) %Get analytic solution at t_{n+1} - fval_analytic = exact_solution_vector(pde,opts,hash_table,t+dt); + + fval_analytic = exact_solution_vector(pde,opts,hash_table,time); %fval_analytic = exact_solution_vector(pde,opts,hash_table,t); %Get analytic moments @@ -728,7 +729,7 @@ subplot(2,2,3); plot(nodes,pde.params.th(nodes)-analytic_moments.T'); title('th_f'); - sgtitle("Fluid Variables. t = "+num2str(t+dt)); + sgtitle("Fluid Variables. t = "+num2str(time)); drawnow else %Plot moments @@ -745,7 +746,7 @@ subplot(2,2,3); plot(nodes,pde.params.th(nodes)); title('th_f'); - sgtitle("Fluid Variables. t = "+num2str(t+dt)); + sgtitle("Fluid Variables. t = "+num2str(time)); drawnow end From cdf6feecc75cf42352d77d6ff38672a67d79f8f2 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Fri, 20 May 2022 08:39:20 -0400 Subject: [PATCH 09/12] commit before branching for limiter --- pdes/vlasov_lb_full_f.m | 6 ++++-- time_advance.m | 33 ++++++++++++++++++++++++++++++--- 2 files changed, 34 insertions(+), 5 deletions(-) diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index 741c5e28..cc9e4e99 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -10,8 +10,7 @@ % % % implicit -% asgard(@vlasov_lb_full_f,'timestep_method','IMEX','deg',3,'lev',[8 3],'dt',0.0002,'num_steps',500,'grid_type','FG','output_grid','interp','quiet',true,'build_realspace_output',false) -% +%asgard(@vlasov_lb_full_f,'timestep_method','IMEX','deg',3,'lev',[6 5],'dt',0.0004,'num_steps',250,'grid_type','SG','output_grid','interp','quiet',false,'build_realspace_output',false,'case',2,'plot_freq',50) %% Define the dimensions % @@ -118,6 +117,9 @@ ic1 = new_md_func(num_dims,{@(x,p,t) (0.5*sin(2*pi*x)+1)./sqrt(2*pi*T_0),... @(v,p,t) exp(-(v-U_0).^2/(2*T_0)),... @(t,p) 0*t+1}); + ic2 = new_md_func(num_dims,{@(x,p,t) 0*x+1,... + @(v,p,t) 0*v+1,... + @(t,p) 0*t+1}); initial_conditions = {ic1}; end diff --git a/time_advance.m b/time_advance.m index 8f4238e4..2bdcf3a0 100644 --- a/time_advance.m +++ b/time_advance.m @@ -512,7 +512,10 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -persistent P +persistent hash_table_FG A_data_FG perm_FG iperm_FG pvec_FG +persistent perm_SG iperm_SG pvec_SG +persistent per iper %Conversion for limiters +persistent FMWT_2D %Switch for IMEX iteration: % 1 : Backward Euler in explicit terms 'E', @@ -529,6 +532,21 @@ %Make fake pde file for use in physical realspace transformation pde_1d.dimensions = pde.dimensions(1); +if isempty(hash_table_FG) + %hash_table_FG = hash_table_nD(pde.get_lev_vec(),'FG'); + [hash_table_FG.elements, hash_table_FG.elements_idx] = hash_table_sparse_nD (pde.get_lev_vec, opts.max_lev, 'FG'); + A_data_FG = global_matrix(pde,opts,hash_table_FG); + [perm_FG,iperm_FG,pvec_FG] = sg_to_fg_mapping_2d(pde,opts,A_data_FG); + [perm_SG,iperm_SG,pvec_SG] = sg_to_fg_mapping_2d(pde,opts,A_data); + x = pde.dimensions{1}.min:(pde.dimensions{1}.max-pde.dimensions{1}.min)/2^pde.dimensions{1}.lev:pde.dimensions{1}.max; + v = pde.dimensions{2}.min:(pde.dimensions{2}.max-pde.dimensions{2}.min)/2^pde.dimensions{2}.lev:pde.dimensions{2}.max; + per = reshape(convertVectoMat(x,v,opts.deg-1,1:(numel(x)-1)*(numel(v)-1)*(opts.deg)^2),[],1); + iper(per) = 1:numel(per); + FMWT_x = OperatorTwoScale_wavelet2(opts.deg,pde.dimensions{1}.lev); + FMWT_v = OperatorTwoScale_wavelet2(opts.deg,pde.dimensions{2}.lev); + FMWT_2D = kron(FMWT_x,FMWT_v); +end + %Create moment matrices that take DG function in (x,v) and transfer it %to DG function in x. @@ -544,7 +562,7 @@ % [~,S,V] = svd(full(M),'econ'); % P = V(:,1:sum(diag(S) > 1e-12)); %end - + hash_table_1D = hash_table_2D_to_1D(hash_table,opts); assert(isempty(pde.termsLHS),'LHS terms currently not supported by IMEX'); @@ -611,8 +629,16 @@ %%%%% %Explicit step - f_2s = f0 + dt*fast_2d_matrix_apply(opts,pde,A_data,f0,'E'); + f_2s = f0 + dt*fast_2d_matrix_apply(opts,pde,A_data,f0,'E'); %f_2s = f0; + %f_2sl = zeros(size(pvec_SG)); f_2sl(pvec_SG) = f_2s(perm_SG(pvec_SG)); f_2sl = f_2sl(iperm_FG); + + %%%%%% Testing slope limiters for SG functions. + %%%%%% C++ do not implement %%%%%%%%% + %f_2sl = limitWrapper(pde,opts,perm_FG,iperm_FG,pvec_FG,perm_SG,iperm_SG,pvec_SG,per,iper,FMWT_2D,f_2s); + %%%%%% END %%%%%% + + %Create rho_2s mom0 = moment_mat{1}*f_2s; %integral of (f,1)_v @@ -664,6 +690,7 @@ f_3s = f0 + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f0+f_2,'E') ... + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f_2,'I'); + f_3sl = limitWrapper(pde,opts,perm_FG,iperm_FG,pvec_FG,perm_SG,iperm_SG,pvec_SG,per,iper,FMWT_2D,f_3s); %f_3s = f0; %Create rho_3s From fdd1d4a91325a4decba2297d755a845990e9d61a Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Fri, 20 May 2022 15:12:26 -0400 Subject: [PATCH 10/12] removing limiting piece from branch --- time_advance.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/time_advance.m b/time_advance.m index 2bdcf3a0..6a83c014 100644 --- a/time_advance.m +++ b/time_advance.m @@ -690,7 +690,7 @@ f_3s = f0 + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f0+f_2,'E') ... + 0.5*dt*fast_2d_matrix_apply(opts,pde,A_data,f_2,'I'); - f_3sl = limitWrapper(pde,opts,perm_FG,iperm_FG,pvec_FG,perm_SG,iperm_SG,pvec_SG,per,iper,FMWT_2D,f_3s); + %f_3sl = limitWrapper(pde,opts,perm_FG,iperm_FG,pvec_FG,perm_SG,iperm_SG,pvec_SG,per,iper,FMWT_2D,f_3s); %f_3s = f0; %Create rho_3s From b0825c00dd7e6673c47185fc38a9523886683ef6 Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Thu, 29 Dec 2022 15:47:14 -0500 Subject: [PATCH 11/12] Added 1x-3v case to moment_reduced_matrix --- moment_reduced_matrix.m | 69 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 62 insertions(+), 7 deletions(-) diff --git a/moment_reduced_matrix.m b/moment_reduced_matrix.m index 9c20a317..8d18ea55 100644 --- a/moment_reduced_matrix.m +++ b/moment_reduced_matrix.m @@ -56,7 +56,7 @@ x_dim = 1; %Hardcoded for now. Needs to change in multi-dim v_dim_1 = 2; %Hardcoded for now. Needs to change in multi-dim - v_dim_2 = 2; %Hardcoded for now. Needs to change in multi-dim + v_dim_2 = 3; %Hardcoded for now. Needs to change in multi-dim %Get row-dim of M -- Assuming a full grid output here rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); @@ -73,21 +73,21 @@ for i=1:num_ele %loop through all elements lit up %Get local dof for element - l_dof = [A_data.element_local_index_D{x_dim }(i)-1,... - A_data.element_local_index_D{v_dim_1}(i)-1,... - A_data.element_local_index_D{v_dim_2}(i)-1]; + l_dof_x = A_data.element_local_index_D{x_dim }(i)-1; + l_dof_v = [A_data.element_local_index_D{v_dim_1}(i)-1,... + A_data.element_local_index_D{v_dim_2}(i)-1]; for xdeg=1:deg %The row of the matrix is given by the x-dof - row_idx = l_dof(1)*deg+xdeg; + row_idx = l_dof_x*deg+xdeg; for vdeg1=1:deg for vdeg2=1:deg %Get column entry based on degree index col_idx = (i-1)*deg^3 + (xdeg-1)*deg^2 + (vdeg1-1)*deg + vdeg2; %Get entry (which is the product of the 1d elements of %g_vec) - val = g_vec_1((A_data.element_local_index_D{v_dim_1}(i)-1)*deg+vdeg1) * ... - g_vec_2((A_data.element_local_index_D{v_dim_2}(i)-1)*deg+vdeg2); + val = g_vec_1(l_dof_v(1)*deg+vdeg1) * ... + g_vec_2(l_dof_v(2)*deg+vdeg2); I(count+1) = row_idx; J(count+1) = col_idx; S(count+1) = val; @@ -99,7 +99,62 @@ %Construct sparse matrix M = sparse(I,J,S,rows,deg^3*num_ele); + + +elseif num_dims == 4 % 1x-3v + + x_dim = 1; %Hardcoded for now. Needs to change in multi-dim + v_dim_1 = 2; %Hardcoded for now. Needs to change in multi-dim + v_dim_2 = 3; %Hardcoded for now. Needs to change in multi-dim + v_dim_3 = 4; %Hardcoded for now. Needs to change in multi-dim + + %Get row-dim of M -- Assuming a full grid output here + rows = deg*2^ceil(log2(max(A_data.element_local_index_D{x_dim}))); + + %Containers for sparse matrix + I = zeros(deg^4*num_ele,1); + J = zeros(deg^4*num_ele,1); + S = zeros(deg^4*num_ele,1); + + g_vec_1 = pde.moments{moment_idx}.fList{v_dim_1}; + g_vec_2 = pde.moments{moment_idx}.fList{v_dim_2}; + g_vec_3 = pde.moments{moment_idx}.fList{v_dim_3}; + + count = 0; + for i=1:num_ele %loop through all elements lit up + + %Get local dof for element + l_dof_x = A_data.element_local_index_D{x_dim }(i)-1; + l_dof_v = [A_data.element_local_index_D{v_dim_1}(i)-1,... + A_data.element_local_index_D{v_dim_2}(i)-1,... + A_data.element_local_index_D{v_dim_3}(i)-1]; + + for xdeg=1:deg + %The row of the matrix is given by the x-dof + row_idx = l_dof_x*deg+xdeg; + for vdeg1=1:deg + for vdeg2=1:deg + for vdeg3=1:deg + %Get column entry based on degree index + col_idx = (i-1)*deg^4 + (xdeg-1)*deg^3 + (vdeg1-1)*deg^2 + (vdeg2-1)*deg + vdeg3; + %Get entry (which is the product of the 1d elements of + %g_vec) + val = g_vec_1(l_dof_v(1)*deg+vdeg1) * ... + g_vec_2(l_dof_v(2)*deg+vdeg2) * ... + g_vec_3(l_dof_v(3)*deg+vdeg3); + I(count+1) = row_idx; + J(count+1) = col_idx; + S(count+1) = val; + count = count + 1; + end + end + end + end + end + %Construct sparse matrix + M = sparse(I,J,S,rows,deg^4*num_ele); + end From 1c9278c6d22ddda42e2986a399abf1dc1698a89f Mon Sep 17 00:00:00 2001 From: Stefan Schnake Date: Thu, 29 Dec 2022 15:48:11 -0500 Subject: [PATCH 12/12] Added to comment --- moment_reduced_matrix.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/moment_reduced_matrix.m b/moment_reduced_matrix.m index 8d18ea55..912bf48d 100644 --- a/moment_reduced_matrix.m +++ b/moment_reduced_matrix.m @@ -2,7 +2,8 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% CURRENTLY ONLY WORKS FOR 1x-1v or 1x-2v system +% CURRENTLY ONLY WORKS FOR 1x-1v, +% 1x-2v, or 1x-3v system %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %