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

[Andrew Spott]

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>