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/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/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..912bf48d 100644 --- a/moment_reduced_matrix.m +++ b/moment_reduced_matrix.m @@ -1,9 +1,10 @@ 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, +% 1x-2v, or 1x-3v system +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % This function produces a matrix M that is the transformation matrix that @@ -15,40 +16,147 @@ % 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 -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 -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'; + 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 = 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}))); + + %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 - temp = g_vec((A_data.element_local_index_D{v_dim}(i)-1)*deg + k); - S(count:count+deg-1) = temp'; + %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]; + + 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 + %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(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; + count = count + 1; + end + end + end + end + + %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 - count = count + deg; + %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 -end + + %Construct sparse matrix + M = sparse(I,J,S,rows,deg^4*num_ele); -M = sparse(I,J,S,rows,deg^2*num_ele); +end end 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 diff --git a/pdes/vlasov_lb_full_f.m b/pdes/vlasov_lb_full_f.m index 35d679a4..cc9e4e99 100644 --- a/pdes/vlasov_lb_full_f.m +++ b/pdes/vlasov_lb_full_f.m @@ -10,12 +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) -% - -soln_x = @(x,p,t) 0*x+1; -soln_v = @(v,p,t) 0*v+1; -soln_t = @(t,p) 0*t+1; +%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 % @@ -52,35 +47,81 @@ %% 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). %soln1 = new_md_func(num_dims,{soln_x,soln_v,soln_t}); %solutions = {soln1}; -solutions = {}; +switch opts.case_ + 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 %% 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) (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 %% Define the terms of the PDE % @@ -152,7 +193,7 @@ % div_v(th q) % q = \grad_v f -g1 = @(x,p,t,dat) 1; +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); 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..6a83c014 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,10 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +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', @@ -526,17 +526,28 @@ 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); + +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. if numel(pde.dimensions) >= 2 @@ -546,6 +557,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 +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 @@ -626,6 +652,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 +682,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 @@ -643,6 +690,8 @@ 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 mom0 = moment_mat{1}*f_3s; %integral of (f,1)_v @@ -666,12 +715,51 @@ 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; - %Plot moments - if ~opts.quiet || 1 + 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,time); + %fval_analytic = exact_solution_vector(pde,opts,hash_table,t); + + %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)-analytic_moments.n'); + title('n_f'); + subplot(2,2,2); + plot(nodes,pde.params.u(nodes)-analytic_moments.u'); + title('u_f'); + subplot(2,2,3); + plot(nodes,pde.params.th(nodes)-analytic_moments.T'); + title('th_f'); + sgtitle("Fluid Variables. t = "+num2str(time)); + drawnow + else + %Plot moments fig1 = figure(1000); fig1.Units = 'Normalized'; @@ -685,11 +773,13 @@ 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 + %fprintf("Norm Check t = %f, |f1|_2 = %f\n",t+dt,norm(f1)); + end end