<div class="gmail_extra">On Tue, May 1, 2012 at 2:27 PM, Andrew Spott <span dir="ltr">&lt;<a href="mailto:andrew.spott@gmail.com" target="_blank">andrew.spott@gmail.com</a>&gt;</span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word">This is actually what I figured, but I haven&#39;t seen the jacobian used in this context before, and wanted to make sure before I bugged you guys with a more complicated problem.<div><br>
</div><div>I am attempting to do a very basic linear time stepping problem:<div><br></div><div>u_t = H(t) u</div><div><br></div><div>Where H(t) = H0 + a(t) * H&#39;</div><div><br></div><div>Where a(t) is just a scalar function of time, and H&#39; has nothing on the diagonal (but is symmetric) and H0 is diagonal</div>
<div><br></div><div>I&#39;m following ex4.c because my problem is very similar to that one, using a linear rhs.</div><div><br></div><div>I create my time step context, and the matrix to use:</div><div><br></div><div>    //creating wave function, creating </div>
<div><br></div><div><div><div>    MatCreate(world, &amp;A);</div><div>    MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,evalues.size(),evalues.size());</div><div>    MatSetType(A,MATMPIAIJ);</div><div>    MatSetFromOptions(A);</div>
</div></div><div><div>    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);     // I don&#39;t know why I need these lines, but when I don&#39;t have them, the program complains about the matrix being in the wrong state.</div><div>
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);</div></div><div><br></div><div><div>    TSCreate(world, &amp;timeStepContext);</div><div>    TSSetProblemType(timeStepContext,TS_LINEAR);</div><div>    TSMonitorSet(timeStepContext,Monitor,&amp;cntx,PETSC_NULL);</div>
<div>    TSSetType(timeStepContext,TSCN);</div><div><br></div><div>    TSSetRHSFunction(timeStepContext,PETSC_NULL,TSComputeRHSFunctionLinear,&amp;cntx);</div><div>    TSSetRHSJacobian(timeStepContext,A,A,*Hamiltonian,&amp;cntx);</div>
<div>    TSSetInitialTimeStep(timeStepContext,0.0,cntx.params-&gt;dt());</div><div>    TSSetSolution(timeStepContext,wavefunction);</div></div><div><br></div><div>Then, my hamiltonian function is:</div><div><br></div><div>
<div>PetscErrorCode Hamiltonian(TS timeStepContext, PetscReal t, Vec u, Mat *A, Mat *B, MatStructure *flag, void* context)</div><div>{</div><div>    Mat h = *A;</div><div>    PetscErrorCode  err;</div><div>    Context         *cntx = (Context*)context;</div>
<div>    PetscScalar     ef = eField(t, cntx-&gt;params);</div><div>    </div><div>    </div><div>    MatDuplicate(cntx-&gt;energy_eigenstates,MAT_COPY_VALUES,&amp;h);</div></div></div></div></blockquote><div><br></div><div>
This is a problem. You are overwriting &#39;h&#39; here, so it is never written back to *A.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">
<div><div><div>    MatAXPY(h, ef, cntx-&gt;dipole_matrix, DIFFERENT_NONZERO_PATTERN);</div><div>    </div><div>    MatAssemblyBegin(h, MAT_FINAL_ASSEMBLY);</div><div>    MatAssemblyEnd(h, MAT_FINAL_ASSEMBLY);</div><div><br>
</div><div>    //If I test the matrix &#39;h&#39; here, against a vector (1,0,0….), I get the correct results that I would expect.</div></div></div></div></blockquote><div><br></div><div>You need</div><div><br></div><div>
     *A = h;</div><div><br></div><div>in order to update it.</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">
<div><div><div>}</div><div> </div></div><div><br></div><div>But this doesn&#39;t work.  my &quot;Monitor&quot; function writes the wave function to disk, but it remains (1,0,0,0…,0,0)</div><div><br></div><div>When I do a test on h inside the Hamiltonian function, (create a vector that is (1,0,0,0…0,0), and multiply it by h), I get the expected results, so the matrices are correct, and h is being created, but when I write out the vector as it is time stepped, (in the monitor function), nothing changes (I get (1,0,0,0…,0,0) for each time step)</div>
<div><br></div><div>My monitor code for reference:</div><div><br></div><div><div>PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec u, void *cntx)</div><div>{</div><div>    Context         *context = (Context*) cntx;</div>
<div>    PetscErrorCode  ierr;</div><div>    PetscScalar     ef = eField(t, context-&gt;params);</div><div>    if (step &lt; 30)</div><div>    {</div><div>        std::stringstream s(&quot;&quot;);</div><div>        s &lt;&lt; step;</div>
<div>        ierr = PetscViewerASCIIOpen(context-&gt;world, context-&gt;params-&gt;getWaveFunctionFilename().append(s.str()).append(&quot;.dat&quot;).c_str(),&amp;context-&gt;view);</div><div>        ierr = PetscViewerSetFormat(context-&gt;view,<span style="white-space:pre-wrap">        </span> PETSC_VIEWER_ASCII_SYMMODU);</div>
<div>        ierr = VecView(u, context-&gt;view);</div><div>    }</div><div>    return ierr;</div><div>}</div></div><div><br></div><div>I&#39;m sorry to bug you guys with this, but I&#39;m a little lost on where even to start looking for an error…</div>
<div><br></div><div>Thanks as always for the help,</div><div><br></div><div>-Andrew</div><div><br></div><div><br></div><div>On May 1, 2012, at 8:36 AM, Matthew Knepley wrote:</div><div><div><br><blockquote type="cite"><div class="gmail_extra">
On Tue, May 1, 2012 at 10:24 AM, Hong Zhang <span dir="ltr">&lt;<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>&gt;</span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

Andrew :<div>We have TS examples under</div><div>~/src/ts/examples/tutorials/</div><div>/src/ts/examples/tests</div><div><br></div><div>Hong<br><div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">


I want to solve a very simple equation:<br>
<br>
u_t = F(t) u<br>
<br>
Where F(t) = H_0 + a(t) H&#39; (H_0 and H&#39; are constant matrices, and a(t) is a time dependent scalar).<br>
<br>
But I&#39;m not sure how to go about doing this using the TS context.<br>
<br>
I don&#39;t have a Jacobian that I need to be worried about, so should I be doing:<br></blockquote></div></div></div></blockquote><div><br></div><div>I am not sure I understand what you mean. Your Jacobian above is F(t), so it does change. You can</div>

<div>of course do this MF since it will be exact, just probably more expensive.</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

<div><div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
TSSetRHSFunction(ts,PETSC_NULL,myRHSFunction,&amp;appctx);<br>
TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&amp;appctx);<br>
<br>
Where:<br>
myRHSFunction(TS ts,PetscReal t,Vec u,Vec F,void *ctx)<br>
{<br>
        //Create temporary matrix A = H_0 + a(t) H&#39;<br>
        //then do F = A u<br>
}<br>
<br>
Or should I be doing something else?<br>
<br>
Thanks for the help, unfortunately, it looks like the documentation on TS in the manual isn&#39;t accurate.<br>
<span><font color="#888888"><br>
-Andrew</font></span></blockquote></div><br></div></div>
</blockquote></div><br><br clear="all"><span class="HOEnZb"><font color="#888888"><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<br>
</font></span></div>
</blockquote></div><br></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<br>
</div>