|
44 | 44 | dz = L/(Nz-1);
|
45 | 45 | r = linspace(0,R,Nr)';
|
46 | 46 | z = linspace(0,L,Nz)';
|
| 47 | +n = Nr*Nz; |
47 | 48 |
|
48 | 49 | % Operational parameters (play around)
|
49 | 50 | vmax = 1;
|
|
54 | 55 | Cf = ones(Nr,1); % Feed concentration
|
55 | 56 |
|
56 | 57 | % The initial condition
|
57 |
| -y0 = zeros(Nr*Nz,1); |
58 |
| -yp0 = model(0,y0,zeros(Nr*Nz,1)); |
| 58 | +y0 = zeros(n,1); |
| 59 | +zz = zeros(n,1); |
| 60 | +yp0 = model(0,y0,zz); |
59 | 61 |
|
60 | 62 | % Auto generate the sparsity pattern (don't try to run without it if Nr > 30 and/or Nz > 30)
|
61 | 63 | % Try to build manually and compare with the auto-generated ;)
|
62 |
| -Jp = zeros(Nr*Nz); |
63 |
| -yr = rand(Nr*Nz,1); |
64 |
| -zz = zeros(Nr*Nz,1); |
| 64 | +% |
| 65 | +tic |
| 66 | +Jp = spalloc(n,n,3*n); |
| 67 | +yr = rand(n,1); |
65 | 68 | delta = 1e-4;
|
66 |
| -for ii = 1:Nr*Nz |
67 |
| - per = zeros(Nr*Nz,1); |
| 69 | +per = zeros(n,1); |
| 70 | +for ii = 1:n |
68 | 71 | per(ii) = 1;
|
69 | 72 | f0 = model(0,yr,zz);
|
70 | 73 | f1 = model(0,yr+per*delta,zz);
|
71 |
| - Jp(:,ii) = (f1 - f0)/delta ~= 0; |
| 74 | + Jp(:,ii) = (f1 - f0)/delta ~= 0; %#ok<SPRIX> |
| 75 | + per(ii) = 0; |
72 | 76 | end
|
73 | 77 |
|
74 |
| -Jp = sparse(Jp); |
75 |
| -Jyp = speye(Nr*Nz); |
| 78 | +Jyp = speye(n); |
| 79 | +etime = toc; |
| 80 | + |
| 81 | +fprintf('Elapsed time building the sparsity pattern = %2.4f s\n',etime); |
| 82 | + |
76 | 83 |
|
77 | 84 | % The solution
|
78 |
| -tf = 25; |
| 85 | +tf = 20; |
79 | 86 | frames = 2000;
|
80 | 87 | tspan = linspace(0,tf,frames);
|
81 | 88 | op = odeset('AbsTol',1e-6,'RelTol',1e-4,'JPattern',{Jp,Jyp});
|
82 | 89 |
|
83 | 90 | tic
|
84 | 91 | [t,y] = ode15i(@model,tspan,y0,yp0,op);
|
85 |
| -toc |
| 92 | +etime = toc; |
| 93 | + |
| 94 | +fprintf('Elapsed time integrating the system = %2.4f s\n',etime); |
| 95 | + |
86 | 96 |
|
87 | 97 | % Parse the solution matrix
|
| 98 | +tic |
88 | 99 | Nt = length(t);
|
89 | 100 | Csol = zeros(Nr,Nz,Nt);
|
90 | 101 | for ii = 1:Nt
|
91 | 102 | Csol(:,:,ii) = reshape(y(ii,:)',Nr,Nz);
|
92 | 103 | end
|
| 104 | +etime = toc; |
| 105 | + |
| 106 | +fprintf('Elapsed time parsing the solution = %2.4f s\n',etime); |
93 | 107 |
|
94 | 108 | % Plot the data (to stop the animation, press CTRL+C on command window)
|
95 | 109 | close all
|
|
137 | 151 | % In z = L and r = 2:Nr-1
|
138 | 152 | res(2:Nr-1,Nz) = C(2:Nr-1,Nz) - C(2:Nr-1,Nz-1);
|
139 | 153 |
|
140 |
| - res = reshape(res,Nr*Nz,1); |
| 154 | + res = reshape(res,n,1); |
141 | 155 |
|
142 | 156 | end
|
143 | 157 |
|
|
0 commit comments