|
On 6/22/2011 1:29 AM, Eranda Gunathilaka wrote:
> "Nasser M. Abbasi"<nma@12000.org> wrote in message<itpbq5$41s$1@speranza.aioe.org>...
>> On 6/20/2011 10:15 PM, Eranda Gunathilaka wrote:
>>
>>> It also seems not fast enough and came to know about the Direct Sparse
>>> Solver for convergence. Can anyone suggest me how to put it in Matlab?
>>>
>>
>> For direct solver, use the slash. It will detect if A is sparse
>> and do the right thing.Matlab I think uses UMFPACK
>> (multifrontal direct solver method) for sparse direct solver.
>>
>> For SOR, it can be slow or fast, depending on how optimal
>> the SOR parameter you choose. How did you choose your SOR
>> parameter? Direct solver will be faster, if you have the
>> needed RAM. Even with sparse matrix, you will run out of
>> ram for larger problems, that is why on real large problems
>> iterative solvers are used much more. For your 500x500, direct
>> solver should not be a problem. This is small problem.
>>
>> --Nasser
>
>
> Hi Nasser,
>
> Thanks for your suggestion.
> Here is my coding for 20x30 grid sample.
> n=20;
> m=30;
> T = zeros(n,m); % All T(i,j) = 0 including the boundaries
> T(1,1)=200;T(15,30)=100; % Known sea Temperature points
> TN = T; % TN = new iteration for solution
> E = TN-T; % Relative error between iterations
> EP = 1e-5; % tolerance for convergence
> imax = 100000; % maximum allowed iterations
> k = 1; % initialization of the iteration
> % Computation
> while k<= imax
> for i = 2:n-1
> for j = 2:m-1
> TN(i,j)=(T(i-1,j)+T(i+1,j)+(T(i,j-1)+T(i,j+1)))/4;% All Center cells
> TN(1,j)=(2*T(2,j)+T(1,j-1)+T(1,j+1))/4;% Bottom Boundary cells
> TN(n,j)=(2*T(n-1,j)+T(n,j-1)+T(n,j+1))/4;% Upper Boundary cells
> TN(i,1)=(2*T(i,2)+T(i-1,1)+T(i+1,1))/4;% Left Boundary cells
> TN(i,m)=(2*T(i,m-1)+T(i-1,m)+T(i+1,m))/4;% Right Boundary cells
> TN(1,1)=(T(1,2)+T(2,1))/2;% Bottom Left Corner Boundary cell
> TN(1,m)=(T(1,m-1)+T(2,m))/2;% Bottom Right Corner Boundary cell
> TN(n,m)=(T(n,m-1)+T(n-1,m))/2;% Upper Right Corner Boundary cell
> TN(n,1)=(T(n-1,1)+T(n,2))/2;% Upper Left Corner Boundary cell
> T(1,1)=200;T(15,30)=100;% Known Temp Values
> E(i,j) = abs(TN(i,j)-T(i,j)); % Relative Error matrix
> end;
> end;
> T = TN;
> k = k + 1;
> Emax = max(max(E));
> if Emax< EP
> figure(1);C=contour(T,20);clabel(C); % Result Plot
> return
> end;
> end;
>
> I used try different values of SOR parameter (w) and 1.85 gave me the lowest iterations. For that I
> modified the statement for all centred cells to:
> .....
> w=1.85;
> .....
> TN(i,j)=((1-w)*T(i,j)+ w*(TN(i-1,j)+TN(i+1,j)+(TN(i,j-1)+T(i,j+1)))/4);
> .....
>
> It reduced the total number of iterations from 17069 to 2580.
> However when I increased the cells to 100x100.
> It is taking very long time.
> My PC is 1G RAM and 2.4 GHz Core 2 DUO Intel Processor.
> I still need to expand to 500x500 for the real case.
> Do you think my PC can handle this.?
>
> Cant I use the same value 1.85 for SOR (w) for this gird as well.?
> I am studying about your proposed method 'Direct Sparse Solver'.
> Do you think the DSS will work for me in my case (still I have no idea how to restructure my coding to DSS).
>
> Looking forward for your valuable suggestions and some materials on DSS for my reading.
>
> Thanks once again.
>
> Eranda.
few general things:
1. To check when to terminate, use relative residual norm. i.e. calculate
the residual R at each iteration along with the solution U,
and check if Norm(Residual)/Norm(f) < epsilon.
f is the initial grid state due to the force (RHS). This is a constant
grid. The residual R = f - (1/h^2)(u(i-1,j)+u(i+1,j)+....-4 u(i,j))
Make sure to use norm grid, not standard norm, i.e. to find the 2-norm,
find the 2-norm of the grid, then multiply by sqrt(h). For 1-norm,
just multiply by h, for max-norm, no change needed.
This is snippet of one program I wrote sometime ago:
---------------------
while not(done)
for i = 2 : nPoints-1
for j = 2 : nPoints-1
resid(i,j)= f(i,j) - 1/h^2 * ( u(i-1,j) + u(i+1,j) + ...
u(i,j-1) + u(i,j+1) - 4*u(i,j) );
unew(i,j) = w/4* ( unew(i-1,j) + u(i+1,j) + unew(i,j-1)+...
u(i,j+1) - h^2*f(i,j)) + (1-w)*u(i,j);
end
end
u = unew;
residv = reshape(resid,nPoints^2,1); % use grid vector norm
normResidue = sqrt(h) * norm(residv,2);
if ( normResidue / normf) <tolerance
done = true;
else
k = k+1;
end
end
------------------------------
where in the above, normf is
-----------------------------
% evaluate f(x,y) on grid
[X,Y] = meshgrid(0:h:1, 0:h:1);
f = -exp(-(X-0.25).^2-(Y-0.6).^2);
nPoints = size(f,1);
fv = reshape(f,nPoints^2,1); % use grid vector norm
normf = sqrt(h) * norm(fv,2);
----------------------------------------
2. It should not take too long to solve 100x100 with SOR. It depends
on your epsilon ofcourse. start with epsilon 0.1 *h^2 where h is the grid
spacing. if it takes more than 1-2 seconds, there is something
wrong.
3. for sparse solver, all what you have to do is make A as sparse.
use
Try this to make A, then solve with the slash.
------------------------
function L2 = lap2d(nx,ny)
Lx=lap1d(nx);
Ly=lap1d(ny);
Ix=speye(nx);
Iy=speye(ny);
L2=kron(Iy,Lx)+kron(Ly,Ix);
end
function L=lap1d(n)
e=ones(n,1);
L=spdiags([e -2*e e],[-1 0 1],n,n);
end
--------------------
I do not play with pde's these days, I am into rigid body
dynamics now, so I am little rusty on more details, but the above
may be will help you.
Here is the complete script to do 2D using SOR from one
of my hw's. use at your own risk. I used an optimal formula
to find omega for sor.
--------------------------------------
% file name: nma_HW3_prob3_parta_SOR.m
%
% This script solves the SOR for 2D for HW3, problem 3, parta
% Math 228A, Fall 2010, UC Davis
% Nasser M. Abbasi
%
close all; clear all;
DOPLOTS = false; %set to false if do not want to see plots
% setup the spacings needed for the problem
spacings = [2^-5 2^-6 2^-7 2^-8 2^-9 2^-10];
omega = arrayfun(@(i) 2*(1-pi*spacings(i)),1:length(spacings));
for m = 1:length(spacings)
h = spacings(m);
% evaluate f(x,y) on grid
[X,Y] = meshgrid(0:h:1, 0:h:1);
f = -exp(-(X-0.25).^2-(Y-0.6).^2);
nPoints = size(f,1);
fv = reshape(f,nPoints^2,1); % use grid vector norm
normf = sqrt(h) * norm(fv,2);
w = omega(m);
fprintf('*******************\n');
fprintf('gridsize=[%d,%d]\n',nPoints-2,nPoints-2);
% initialize space (grid) for residual calculation and for solution
resid = zeros(nPoints,nPoints);
u = zeros(nPoints,nPoints);
unew = u;
done = false; %flag set to true in loop below when it converges
tolerance = h^2; % set tolerance
k = 1; % initialize iteration counter
t = cputime;
while not(done)
for i = 2 : nPoints-1
for j = 2 : nPoints-1
resid(i,j)= f(i,j) - 1/h^2 * ( u(i-1,j) + u(i+1,j) + ...
u(i,j-1) + u(i,j+1) - 4*u(i,j) );
unew(i,j) = w/4* ( unew(i-1,j) + u(i+1,j) + unew(i,j-1)+...
u(i,j+1) - h^2*f(i,j)) + (1-w)*u(i,j);
end
end
u = unew;
residv = reshape(resid,nPoints^2,1); % use grid vector norm
normResidue = sqrt(h) * norm(residv,2);
if ( normResidue / normf) <tolerance
done = true;
else
k = k+1;
end
if DOPLOTS
subplot(1,2,1);
mesh(X ,Y ,resid );
subplot(1,2,2);
mesh(X ,Y ,u );
drawnow;
end
end
fprintf('cpu time for SOR=%f seconds\n',cputime-t);
fprintf('number of iterations = %d\n',k);
end
-------------------------------------
--Nasser
|
Feed for this Thread