<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, Aug 6, 2014 at 6:50 PM, Francis Poulin <span dir="ltr"><<a href="mailto:fpoulin@uwaterloo.ca" target="_blank">fpoulin@uwaterloo.ca</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>
<div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">Hello,<br>
<br>
I am trying to figure out a faster way of doing in petsc4py and thought I would ask to see how this should be done in petsc. If I can figure that out then I think I can figure out how to translate it but I am having issues with the basic algorithm I should
be using.<br>
<br>
I am solving a linear system of ODEs of the form <br>
<br>
B du/dt = A u<br>
<div><br>
<div style="font-family:Tahoma;font-size:13px">
<div style="font-family:Tahoma;font-size:13px">using Crank-Nicholson. So the update equation can be written as<br>
<br>
(B - dt/2*A)*u^(n+1) = (B + dt/2*A)*u^n<br>
<br>
or, if you don't mind using the backslash operator,<br>
<br>
u^(n+1) = (B - dt/2*A)\(B + dt/2*A)*u^n<br>
<br>
The integration is faster of course after I build the matrix <br>
<br>
C = (B - dt/2*A)\(B + dt/2*A)<br>
<br>
however building this is very slow and I am not sure how to build it efficiently. At the moment the majority of my code is just doing that.<br></div></div></div></div></div></blockquote><div><br></div><div>1) For this to make sense, you have to be doing many many more than N timesteps, where N is the number of unknowns. Why</div>
<div> not just solve this equation at each timestep?</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">
<div><div style="font-family:Tahoma;font-size:13px"><div style="font-family:Tahoma;font-size:13px">
In petsc4py I do the following: <br>
<br>
(you can find the full code at <a href="https://github.com/francispoulin/qg_1L_ts/blob/master/qg_1L_ts.py" target="_blank">https://github.com/francispoulin/qg_1L_ts/blob/master/qg_1L_ts.py</a>)<br>
<br>
<div><span> AeT</span> <span>
=</span> <span>PETSc</span><span>.</span><span>Mat</span><span>()</span><span>.</span><span>createAIJ</span><span>([</span><span>Ny</span><span>-</span><span>1</span><span>,</span><span>Ny</span><span>-</span><span>1</span><span>])</span></div>
</div></div></div></div></div></blockquote><div><br></div><div>2) Storing the inverse in a sparse matrix does not make sense. Its likely dense. Also, you have not</div><div> preallocated so it will be slow. Use a dense matrix and both problems will go away.</div>
<div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div style="direction:ltr;font-family:Tahoma;color:#000000;font-size:10pt">
<div><div style="font-family:Tahoma;font-size:13px"><div style="font-family:Tahoma;font-size:13px">
<div> <span>AeT</span><span>.</span><span>setUp</span><span>();</span>
<span>AeT</span><span>.</span><span>assemble</span><span>()</span></div>
<div> <span>ksp</span> <span>
=</span> <span>PETSc</span><span>.</span><span>KSP</span><span>()</span><span>.</span><span>create</span><span>(</span><span>PETSc</span><span>.</span><span>COMM_WORLD</span><span>)</span></div>
<div> <span>ksp</span><span>.</span><span>setOperators</span><span>(</span><span>B</span><span>-</span><span>dto2</span><span>*</span><span>A</span><span>)</span></div>
<div> <span>ksp</span><span>.</span><span>setTolerances</span><span>(</span><span>1e-16</span><span>)</span></div>
<div> <span>pc</span> <span>
=</span> <span>ksp</span><span>.</span><span>getPC</span><span>()</span></div>
<div> <span>pc</span><span>.</span><span>setType</span><span>(</span><span>'none'</span><span>)</span></div>
<div> <span>ksp</span><span>.</span><span>setFromOptions</span><span>()</span></div>
<div><br>
</div>
<div> <span>btemp</span> <span>
=</span> <span>B</span><span>+</span><span>dto2</span><span>*</span><span>A</span></div>
<div> <span>btemp</span><span>.</span><span>assemble</span><span>()</span></div>
<div><br>
</div>
<div> <span>bcol</span> <span>
=</span> <span>PETSc</span><span>.</span><span>Vec</span><span>()</span><span>.</span><span>createMPI</span><span>(</span><span>Ny</span><span>-</span><span>1</span><span>)</span></div>
<div> <span>xcol</span> <span>
=</span> <span>PETSc</span><span>.</span><span>Vec</span><span>()</span><span>.</span><span>createMPI</span><span>(</span><span>Ny</span><span>-</span><span>1</span><span>)</span></div>
<div> <span>bcol</span><span>.</span><span>setUp</span><span>();</span>
<span>xcol</span><span>.</span><span>setUp</span><span>()</span></div>
<div><br>
</div>
<div> <span>sc</span><span>,</span><span>ec</span>
<span>=</span> <span>bcol</span><span>.</span><span>getOwnershipRange</span><span>()</span></div>
<div> <span>for</span> <span>
i</span> <span>in</span> <span>range</span><span>(</span><span>0</span><span>,</span><span>Ny</span><span>-</span><span>1</span><span>):</span></div>
<div> <span>bcol</span><span>[</span><span>sc</span><span>:</span><span>ec</span><span>]</span>
<span>=</span> <span>btemp</span><span>[</span><span>sc</span><span>:</span><span>ec</span><span>,</span><span>i</span><span>]</span></div>
<div> <span>bcol</span><span>.</span><span>assemble</span><span>();</span>
<span>xcol</span><span>.</span><span>assemble</span><span>()</span></div>
<div> <span>ksp</span><span>.</span><span>solve</span><span>(</span><span>bcol</span><span>,</span><span>xcol</span><span>)</span></div>
<div> <span>AeT</span><span>[</span><span>sc</span><span>:</span><span>ec</span><span>,</span><span>i</span><span>]</span>
<span>=</span> <span>xcol</span><span>[</span><span>sc</span><span>:</span><span>ec</span><span>]</span></div>
<div><br>
</div>
<div> <span>AeT</span><span>.</span><span>assemble</span><span>()</span></div>
<br>
We are building the matrix one row at a time and, as accurate as that maybe, it is very slow. Is there perhaps a fast way of doing this kind of operation in petsc?<br>
<br>
Sorry for the bother but any advice would be greatly appreciated.<br>
<br>
Cheers, Francis<br>
<br>
------------------<br>
Francis Poulin <br>
Associate Professor<br>
Associate Chair, Undergraduate Studies <br>
Department of Applied Mathematics<br>
University of Waterloo<br>
<br>
email: <a href="mailto:fpoulin@uwaterloo.ca" target="_blank">fpoulin@uwaterloo.ca</a><br>
Web: <a href="https://uwaterloo.ca/poulin-research-group/" target="_blank">https://uwaterloo.ca/poulin-research-group/</a><br>
Telephone: <a href="tel:%2B1%20519%20888%204567%20x32637" value="+15198884567" target="_blank">+1 519 888 4567 x32637</a><br>
<br>
</div>
</div>
</div>
</div>
</div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>