diff --git a/asgard_run_pde.m b/asgard_run_pde.m index fdb36ca6..d0592ed5 100644 --- a/asgard_run_pde.m +++ b/asgard_run_pde.m @@ -536,7 +536,7 @@ else figs.mass = figure('Name','mass(t)','Units','normalized','Position',[0.7,0.5,0.2,0.3]); end - plot(outputs.mass_t/outputs.mass_t(1)); + plot(mass_t/mass_t(1)); xtitle='timestep'; ytitle='mass/mass(t=0)'; ylim([0,2]); diff --git a/sg_to_fg_mapping_2d.m b/sg_to_fg_mapping_2d.m new file mode 100644 index 00000000..47935d74 --- /dev/null +++ b/sg_to_fg_mapping_2d.m @@ -0,0 +1,76 @@ +function [perm,iperm,pvec] = sg_to_fg_mapping_2d(pde,opts,A_data) +% Maps the sparse_grid index to standard tensor product full_grid index. +% This to allow multiplication by the full-grid kronecker sum stiffness +% matrix. + +% perm -> SG-type index to FG-standard index +% pvec -> logical vector showing which elements are lit. Used for SG->FG +% conversion. +% iperm -> inverse mapping + +% This reindexing allows us to use the standard sum of kroncker products +% stiffness matrix in lieu of the construction used in apply_A. Using the +% conversion to full-grid and back is also much faster. + +% Works for adaptiviity + +% Only works in 2D for now, but can be easily changed for higher diminsions + +% Function assumes A_data.element_global_row_index is 1:sparse_dof. +% --Does not check for this + +% USAGE: If the matrix A and vector f are in the SG-index, then +% we can compute A*f=f1 using A_F by the sequence +% f_F = zeros(size(A_F,1),1); +% f_F(pvec) = f(perm(pvec)); +% Af_F = A_F*f_F; +% f1 = Af_F(iperm); +% where A_F is the standard full-grid sum of kronecker products. +% Specifically, A_F is generated by the code: +% A_F = 0; +% for i=1:numel(pde.terms) +% A_F = A_F + kron(pde.terms{i}.terms_1D{1}.mat,... +% pde.terms{i}.terms_1D{2}.mat); +% end +% +% The SG-index stiffness matrix and be computed from the full grid matrix +% by using: +% A = A_F(iperm,iperm); + +%Checking for only 2 dimensions +assert(numel(pde.dimensions) == 2); + +deg = opts.deg; + +num_sparse_ele = numel(A_data.element_local_index_D{1}); %Number of sparse elements +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); +perm = nan(max_x_dof*max_y_dof,1); + +idx = 1; %Keep track of index +for i=1:num_sparse_ele + %First compute the x dim FG dof + dof_x = (A_data.element_local_index_D{1}(i)-1)*deg + (1:deg); + %Now y dim FG dof + dof_y = (A_data.element_local_index_D{2}(i)-1)*deg + (1:deg); + + %Mapping is (i,j) -> (i-1)*max_y_dof+j; + + for deg_x=1:deg + for deg_y=1:deg + iperm(idx) = (dof_x(deg_x)-1)*max_y_dof+dof_y(deg_y); + idx = idx + 1; + end + end +end + + +perm(iperm) = 1:numel(iperm); +pvec = ~isnan(perm); + + +end + diff --git a/time_advance.m b/time_advance.m index df332238..0db34abb 100644 --- a/time_advance.m +++ b/time_advance.m @@ -251,162 +251,28 @@ s0 = source_vector(pde,opts,hash_table,t+dt); bc0 = boundary_condition_vector(pde,opts,hash_table,t+dt); -if opts.time_independent_A || opts.time_independent_build_A - persistent A; +%Get SG <-> FG conversion +[~,iperm,~] = sg_to_fg_mapping_2d(pde,opts,A_data); + +%Construct full-grid A and ALHS +A_F = 0; +for i=1:numel(pde.terms) + A_F = A_F + kron(pde.terms{i}.terms_1D{1}.mat,... + pde.terms{i}.terms_1D{2}.mat); end - -if opts.time_independent_A % relies on inv() so is no good for poorly conditioned problems - persistent AA_inv; - persistent ALHS_inv; - - if isempty(AA_inv) - if applyLHS - [~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg); - I = speye(numel(diag(A))); - ALHS_inv = inv(ALHS); - AA = I - dt*(ALHS_inv * A); - b = f0 + dt*(ALHS_inv * (s0 + bc0)); - else - [~,A] = apply_A(pde,opts,A_data,f0,deg); - I = speye(numel(diag(A))); - AA = I - dt*A; - b = f0 + dt*(s0 + bc0); - end - - if numel(AA(:,1)) <= 4096 - condAA = condest(AA); - if ~opts.quiet; disp([' condest(AA) : ', num2str(condAA,'%.1e')]); end - - if condAA > 1e6 - disp(['WARNING: Using time_independent_A=true for poorly conditioned system not recommended']); - disp(['WARNING: cond(A) = ', num2str(condAA,'%.1e')]); - end - end - - AA_inv = inv(AA); - f1 = AA_inv * b; - else - if applyLHS - b = f0 + dt*(ALHS_inv * (s0 + bc0)); - else - b = f0 + dt*(s0 + bc0); - end - f1 = AA_inv * b; - end - -else % use the backslash operator instead - if applyLHS - if opts.time_independent_build_A - if isempty(A);[~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg);end - else - [~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg); - end - I = speye(numel(diag(A))); - AA = I - dt*(ALHS \ A); - b = f0 + dt*(ALHS \ (s0 + bc0)); - else - if opts.time_independent_build_A - if isempty(A);[~,A] = apply_A(pde,opts,A_data,f0,deg);end - else - [~,A] = apply_A(pde,opts,A_data,f0,deg); - end - I = speye(numel(diag(A))); - AA = I - dt*A; - b = f0 + dt*(s0 + bc0); - end - - % Direct solve - rescale = false; - if rescale - % [AA_rescaled,diag_scaling] = rescale2(AA); - % AA_thresholded = sparsify(AA_rescaled,1e-5); - [P,R,C] = equilibrate(AA); - AA_rescaled = R*P*AA*C; - b_rescaled = R*P*b; - f1 = AA_rescaled \ b_rescaled; - f1 = C*f1; - else - f1 = AA \ b; +A = A_F(iperm,iperm); +ALHS = speye(size(A)); +if applyLHS + ALHS_F = 0; + for i=1:numel(pde.termsLHS) + ALHS_F = ALHS_F + kron(pde.termsLHS{i}.terms_1D{1}.mat,... + pde.terms{i}.terms_1D{2}.mat); end - - % % Pre-kron solve - % num_terms = numel(pde.terms); - % num_dims = numel(pde.dimensions); - % for d=1:num_dims - % dim_mat_list{d} = zeros(size(pde.terms{1}.terms_1D{1}.mat)); - % for t=1:num_terms - % dim_mat_list{d} = dim_mat_list{d} + pde.terms{t}.terms_1D{d}.mat; - % end - % end - % for d=1:num_dims - % A = dim_mat_list{d}; - % I = speye(numel(diag(A))); - % AA_d = I - dt*A; - % dim_mat_inv_list{d} = inv(AA_d); - % end - % use_kronmultd = true; - % if use_kronmultd - % f1a = kron_multd(num_dims,dim_mat_inv_list,b); - % else - % f1a = kron_multd_full(num_dims,dim_mat_inv_list,b); - % end - - % % Iterative solve - % restart = []; - % tol=1e-6; - % maxit=1000; - % tic; - % [f10,flag,relres,iter,resvec] = gmres(AA,b,restart,tol,maxit); - % t_gmres = toc; - % figure(67) - % semilogy(resvec); - % - % % Direct solve - LU approach to reducing repeated work - % [L,U,P] = lu(AA); - % n=size(AA,1); - % ip = P*reshape(1:n,n,1); % generate permutation vector - % err = norm(AA(ip,:)-L*U,1); % err should be small - % tol = 1e-9; - % isok = (err <= tol * norm(AA,1) ); % just a check - % disp(['isok for LU: ', num2str(isok)]); - % % to solve a linear system A * x = b, we have P * A * x = P*b - % % then from LU factorization, we have (L * U) * x = (P*b), so x = U \ (L \ ( P*b)) - % f1_LU = U \ (L \ (P*b)); - % - % % Direct solve - QR reduced rank - % n=size(AA,1); - % [Q,R,P] = qr(AA); % A*P = Q*R - % % where R is upper triangular, Q is orthogonal, Q’*Q is identity, P is column permutation - % err = norm( AA*P - Q*R,1); - % tol=1e-9; - % isok = (err <= tol * norm(AA,1)); % just a check - % disp(['isok for QR: ', num2str(isok)]); - % % to solve A * x = b, we have A * P * (P’*x) = b, (Q*R) * (P’*x) = b - % % y = R\(Q’*b), P*y = x - % tol=1e-9; - % is_illcond = abs(R(n,n)) <= tol * abs(R(1,1)); - % if(is_illcond) - % disp('is_illcond == true'); - % bhat = Q'*b; - % y = zeros(n,1); - % yhat = R(1:(n-1), 1:(n-1)) \ bhat(1:(n-1)); % solve with smaller system - % y(1:(n-1)) = yhat(1:(n-1)); - % y(n) = 0; % force last component to be zero - % f1_QR = P * y; - % else - % disp('is_illcond == false'); - % f1_QR = P * (R \ (Q'*b)); - % end - % - % disp(['f1-f10: ',num2str(norm(f1-f10)/norm(f1))]); - % disp(['f1-f1_LU: ',num2str(norm(f1-f1_LU)/norm(f1))]); - % disp(['f1-f1_QR: ',num2str(norm(f1-f1_QR)/norm(f1))]); - % disp(['direct runtime: ', num2str(t_direct)]); - % disp(['gmres runtime: ', num2str(t_gmres)]); - % - % f1 = f1_QR; - + ALHS = ALHS_F(iperm,iperm); end + +f1 = (ALHS-dt*A)\(ALHS*f0 + dt*(s0+bc0)); + end