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.
I am attempting to do a very basic linear time stepping problem:
u_t = H(t) u
Where H(t) = H0 + a(t) * H'
Where a(t) is just a scalar function of time, and H' has nothing on the
diagonal (but is symmetric) and H0 is diagonal
I'm following ex4.c because my problem is very similar to that one, using a
linear rhs.
I create my time step context, and the matrix to use:
//creating wave function, creating
MatCreate(world, &A);
MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,evalues.size(),evalues.size());
MatSetType(A,MATMPIAIJ);
MatSetFromOptions(A);
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.
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
TSCreate(world, &timeStepContext);
TSSetProblemType(timeStepContext,TS_LINEAR);
TSMonitorSet(timeStepContext,Monitor,&cntx,PETSC_NULL);
TSSetType(timeStepContext,TSCN);
TSSetRHSFunction(timeStepContext,PETSC_NULL,TSComputeRHSFunctionLinear,&cntx);
TSSetRHSJacobian(timeStepContext,A,A,*Hamiltonian,&cntx);
TSSetInitialTimeStep(timeStepContext,0.0,cntx.params->dt());
TSSetSolution(timeStepContext,wavefunction);
Then, my hamiltonian function is:
PetscErrorCode Hamiltonian(TS timeStepContext, PetscReal t, Vec u, Mat *A, Mat
*B, MatStructure *flag, void* context)
{
Mat h = *A;
PetscErrorCode err;
Context *cntx = (Context*)context;
PetscScalar ef = eField(t, cntx->params);
MatDuplicate(cntx->energy_eigenstates,MAT_COPY_VALUES,&h);
MatAXPY(h, ef, cntx->dipole_matrix, DIFFERENT_NONZERO_PATTERN);
MatAssemblyBegin(h, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(h, MAT_FINAL_ASSEMBLY);
//If I test the matrix 'h' here, against a vector (1,0,0?.), I get the
correct results that I would expect.
}
But this doesn't work. my "Monitor" function writes the wave function to disk,
but it remains (1,0,0,0?,0,0)
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)
My monitor code for reference:
PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec u, void *cntx)
{
Context *context = (Context*) cntx;
PetscErrorCode ierr;
PetscScalar ef = eField(t, context->params);
if (step < 30)
{
std::stringstream s("");
s << step;
ierr = PetscViewerASCIIOpen(context->world,
context->params->getWaveFunctionFilename().append(s.str()).append(".dat").c_str(),&context->view);
ierr = PetscViewerSetFormat(context->view,
PETSC_VIEWER_ASCII_SYMMODU);
ierr = VecView(u, context->view);
}
return ierr;
}
I'm sorry to bug you guys with this, but I'm a little lost on where even to
start looking for an error?
Thanks as always for the help,
-Andrew
On May 1, 2012, at 8:36 AM, Matthew Knepley wrote:
> On Tue, May 1, 2012 at 10:24 AM, Hong Zhang <hzhang at mcs.anl.gov> wrote:
> Andrew :
> We have TS examples under
> ~/src/ts/examples/tutorials/
> /src/ts/examples/tests
>
> Hong
> I want to solve a very simple equation:
>
> u_t = F(t) u
>
> Where F(t) = H_0 + a(t) H' (H_0 and H' are constant matrices, and a(t) is a
> time dependent scalar).
>
> But I'm not sure how to go about doing this using the TS context.
>
> I don't have a Jacobian that I need to be worried about, so should I be doing:
>
> I am not sure I understand what you mean. Your Jacobian above is F(t), so it
> does change. You can
> of course do this MF since it will be exact, just probably more expensive.
>
> Matt
>
> TSSetRHSFunction(ts,PETSC_NULL,myRHSFunction,&appctx);
> TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&appctx);
>
> Where:
> myRHSFunction(TS ts,PetscReal t,Vec u,Vec F,void *ctx)
> {
> //Create temporary matrix A = H_0 + a(t) H'
> //then do F = A u
> }
>
> Or should I be doing something else?
>
> Thanks for the help, unfortunately, it looks like the documentation on TS in
> the manual isn't accurate.
>
> -Andrew
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120501/f4a4adca/attachment.htm>
