[petsc-users] Solving a very simple time step problem:

[Matthew Knepley]

On Tue, May 1, 2012 at 2:27 PM, Andrew Spott <andrew.spott at gmail.com> wrote:

> 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);
>

This is a problem. You are overwriting 'h' here, so it is never written
back to *A.


>     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.
>

You need

     *A = h;

in order to update it.

   Matt


> }
>
>
> 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
>
>
>


-- 
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/b630f4cf/attachment.htm>