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

[Andrew Spott]

Thanks!

-Andrew


On May 1, 2012, at 2:12 PM, Matthew Knepley wrote:

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