<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">This is actually what I figured, but I haven'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'</div><div><br></div><div>Where a(t) is just a scalar function of time, and H' has nothing on the diagonal (but is symmetric) and H0 is diagonal</div><div><br></div><div>I'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>&nbsp; &nbsp; //creating wave function, creating&nbsp;</div><div><br></div><div><div><div>&nbsp; &nbsp; MatCreate(world, &amp;A);</div><div>&nbsp; &nbsp; MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,evalues.size(),evalues.size());</div><div>&nbsp; &nbsp; MatSetType(A,MATMPIAIJ);</div><div>&nbsp; &nbsp; MatSetFromOptions(A);</div></div></div><div><div>&nbsp; &nbsp; MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); &nbsp; &nbsp; // I don't know why I need these lines, but when I don't have them, the program complains about the matrix being in the wrong state.</div><div>&nbsp; &nbsp; MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);</div></div><div><br></div><div><div>&nbsp; &nbsp; TSCreate(world, &amp;timeStepContext);</div><div>&nbsp; &nbsp; TSSetProblemType(timeStepContext,TS_LINEAR);</div><div>&nbsp; &nbsp; TSMonitorSet(timeStepContext,Monitor,&amp;cntx,PETSC_NULL);</div><div>&nbsp; &nbsp; TSSetType(timeStepContext,TSCN);</div><div><br></div><div>&nbsp; &nbsp; TSSetRHSFunction(timeStepContext,PETSC_NULL,TSComputeRHSFunctionLinear,&amp;cntx);</div><div>&nbsp; &nbsp; TSSetRHSJacobian(timeStepContext,A,A,*Hamiltonian,&amp;cntx);</div><div>&nbsp; &nbsp; TSSetInitialTimeStep(timeStepContext,0.0,cntx.params-&gt;dt());</div><div>&nbsp; &nbsp; 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>&nbsp; &nbsp; Mat h = *A;</div><div>&nbsp; &nbsp; PetscErrorCode &nbsp;err;</div><div>&nbsp; &nbsp; Context &nbsp; &nbsp; &nbsp; &nbsp; *cntx = (Context*)context;</div><div>&nbsp; &nbsp; PetscScalar &nbsp; &nbsp; ef = eField(t, cntx-&gt;params);</div><div>&nbsp;&nbsp; &nbsp;</div><div>&nbsp; &nbsp;&nbsp;</div><div>&nbsp; &nbsp; MatDuplicate(cntx-&gt;energy_eigenstates,MAT_COPY_VALUES,&amp;h);</div><div>&nbsp; &nbsp; MatAXPY(h, ef, cntx-&gt;dipole_matrix, DIFFERENT_NONZERO_PATTERN);</div><div>&nbsp; &nbsp;&nbsp;</div><div>&nbsp; &nbsp; MatAssemblyBegin(h, MAT_FINAL_ASSEMBLY);</div><div>&nbsp; &nbsp; MatAssemblyEnd(h, MAT_FINAL_ASSEMBLY);</div><div><br></div><div>&nbsp; &nbsp; //If I test the matrix 'h' here, against a vector (1,0,0….), I get the correct results that I would expect.</div><div>}</div><div>&nbsp;</div></div><div><br></div><div>But this doesn't work. &nbsp;my "Monitor" 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>&nbsp; &nbsp; Context &nbsp; &nbsp; &nbsp; &nbsp; *context = (Context*) cntx;</div><div>&nbsp; &nbsp; PetscErrorCode &nbsp;ierr;</div><div>&nbsp; &nbsp; PetscScalar &nbsp; &nbsp; ef = eField(t, context-&gt;params);</div><div>&nbsp;&nbsp; &nbsp;if (step &lt; 30)</div><div>&nbsp; &nbsp; {</div><div>&nbsp; &nbsp; &nbsp; &nbsp; std::stringstream s("");</div><div>&nbsp; &nbsp; &nbsp; &nbsp; s &lt;&lt; step;</div><div>&nbsp; &nbsp; &nbsp; &nbsp; ierr = PetscViewerASCIIOpen(context-&gt;world, context-&gt;params-&gt;getWaveFunctionFilename().append(s.str()).append(".dat").c_str(),&amp;context-&gt;view);</div><div>&nbsp; &nbsp; &nbsp; &nbsp; ierr = PetscViewerSetFormat(context-&gt;view,<span class="Apple-tab-span" style="white-space: pre; ">        </span>&nbsp;PETSC_VIEWER_ASCII_SYMMODU);</div><div>&nbsp; &nbsp; &nbsp; &nbsp; ierr = VecView(u, context-&gt;view);</div><div>&nbsp; &nbsp; }</div><div>&nbsp; &nbsp; return ierr;</div><div>}</div></div><div><br></div><div>I'm sorry to bug you guys with this, but I'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 class="Apple-interchange-newline"><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&nbsp;:<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' (H_0 and H' are constant matrices, and a(t) is a time dependent scalar).<br>
<br>
But I'm not sure how to go about doing this using the TS context.<br>
<br>
I don'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>&nbsp; &nbsp;Matt</div><div>&nbsp;</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>
 &nbsp; &nbsp; &nbsp; &nbsp;//Create temporary matrix A = H_0 + a(t) H'<br>
 &nbsp; &nbsp; &nbsp; &nbsp;//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't accurate.<br>
<span><font color="#888888"><br>
-Andrew</font></span></blockquote></div><br></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>
</blockquote></div><br></div>
</div></body></html>