Macro micro

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));

asgard_run_pde.m

@@ -89,50 +89,46 @@
 
 %% 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
 
     %%
     % 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;

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

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,36 +17,46 @@
 %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;
 %     end
 % 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