% BVP_field_Jpattern.m % This MATLAB program demonstrates the sparsity patterns % obtained for a 1-D transport-reaction problem % based on the use of different stacking schemes. % K.J. Beers. MIT ChE. 11/15/2004 function iflag = BVP_field_Jpattern(N_f,N_g); iflag = 0; N_tot = N_f*N_g; % total # of unknowns Nz_per_row = N_f + 2; % # of non-zero diagonals Nz = N_tot*Nz_per_row; % # of non-zero elements % scheme I, stacking c1(z1),...,c1(zN),c2(z1),... S_1 = spalloc(N_tot,N_tot,Nz); for i_g = 1:N_g for j_f = 1:N_f m = get_label_1(j_f,i_g,N_f,N_g); % each species at same point for k_f = 1:N_f n = get_label_1(k_f,i_g,N_f,N_g); S_1(m,n) = 1; end % same field, left point if(i_g > 1) n = get_label_1(j_f,i_g-1,N_f,N_g); S_1(m,n) = 1; end % same field, right point if(i_g < N_g) n = get_label_1(j_f,i_g+1,N_f,N_g); S_1(m,n) = 1; end end end % make spy plot of sparsity pattern of Jacobian figure; subplot(1,2,1); spy(S_1); title('Scheme I'); % scheme II, stacking c1(z1),c2(z1),...,c1(z2),c2(z2),... S_2 = spalloc(N_tot,N_tot,Nz); for i_g=1:N_g for j_f = 1:N_f m = get_label_2(j_f,i_g,N_f,N_g); % each species at same point for k_f = 1:N_f n = get_label_2(k_f,i_g,N_f,N_g); S_2(m,n) = 1; end % same field, left point if(i_g > 1) n = get_label_2(j_f,i_g-1,N_f,N_g); S_2(m,n) = 1; end % same field, right point if(i_g < N_g) n = get_label_2(j_f,i_g+1,N_f,N_g); S_2(m,n) = 1; end end end % make spy plot of sparsity pattern of Jacobian subplot(1,2,2); spy(S_2); title('Scheme II'); iflag = 1; return; % ======================== function label = get_label_1(j_f,i_g,N_f,N_g); label = (j_f-1)*N_g + i_g; return; % ======================== function label = get_label_2(j_f,i_g,N_f,N_g); label = (i_g-1)*N_g + j_f; return;