<div class="gmail_extra">On Tue, May 1, 2012 at 2:27 PM, Andrew Spott <span dir="ltr"><<a href="mailto:andrew.spott@gmail.com" target="_blank">andrew.spott@gmail.com</a>></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'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> //creating wave function, creating </div>
<div><br></div><div><div><div> MatCreate(world, &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'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>
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);</div></div><div><br></div><div><div> TSCreate(world, &timeStepContext);</div><div> TSSetProblemType(timeStepContext,TS_LINEAR);</div><div> TSMonitorSet(timeStepContext,Monitor,&cntx,PETSC_NULL);</div>
<div> TSSetType(timeStepContext,TSCN);</div><div><br></div><div> TSSetRHSFunction(timeStepContext,PETSC_NULL,TSComputeRHSFunctionLinear,&cntx);</div><div> TSSetRHSJacobian(timeStepContext,A,A,*Hamiltonian,&cntx);</div>
<div> TSSetInitialTimeStep(timeStepContext,0.0,cntx.params->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->params);</div><div> </div><div> </div><div> MatDuplicate(cntx->energy_eigenstates,MAT_COPY_VALUES,&h);</div></div></div></div></blockquote><div><br></div><div>
This is a problem. You are overwriting 'h' 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->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 'h' 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't work. 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> Context *context = (Context*) cntx;</div>
<div> PetscErrorCode ierr;</div><div> PetscScalar ef = eField(t, context->params);</div><div> if (step < 30)</div><div> {</div><div> std::stringstream s("");</div><div> s << step;</div>
<div> ierr = PetscViewerASCIIOpen(context->world, context->params->getWaveFunctionFilename().append(s.str()).append(".dat").c_str(),&context->view);</div><div> ierr = PetscViewerSetFormat(context->view,<span style="white-space:pre-wrap"> </span> PETSC_VIEWER_ASCII_SYMMODU);</div>
<div> ierr = VecView(u, context->view);</div><div> }</div><div> 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><blockquote type="cite"><div class="gmail_extra">
On Tue, May 1, 2012 at 10:24 AM, Hong Zhang <span dir="ltr"><<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>></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' (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> 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,&appctx);<br>
TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&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'<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'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>