[petsc-users] Using MatSetValuesBlocked to build Finite Element Matrices with variable coefficients

[Nidish]

Update: I had made some mistakes in the assembly in the previous 
version. They're fixed now (see below).

On 8/9/20 9:13 PM, Nidish wrote:
>
> I'm trying to write a simple 1D finite element test case (linear 
> elasticity). It was given in the manual that using an appropriate DM 
> and calling MatSetValuesBlocked(). I'm however unable to do this for 
> cases where I want to take runtime inputs for physical constants in 
> the matrices.
>
> One thing that puzzles me is the use of "const PetscInt*" and "const 
> PetscScalar*" in the declaration of MatSetValuesBlocked(). I can't 
> think of why this is done (if I am to loop through the elements and 
> call this function, I don't see why the indices and/or values must be 
> constant). The declaration I found was:
>
>     PetscErrorCode  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode>  MatSetValues  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/MatSetValues.html#MatSetValues>(Mat  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/Mat.html#Mat>  mat,PetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  m,constPetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  idxm[],PetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  n,constPetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  idxn[],constPetscScalar  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscScalar.html#PetscScalar>  v[],InsertMode  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/InsertMode.html#InsertMode>  addv)
>
> Here's the DM setup for my code:
>
>     DM mshdm;
>     Mater steel = {1.0, 2e11, 7800, 3.0, 400.0};// A structure with
>     {Area, Modulus, Density, Length, Force} (runtime vars)
>     KSP ksp;
>
>     DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL,
>     &mshdm);// N is a runtime variable
>     DMSetFromOptions(mshdm);
>     DMSetUp(mshdm);
>     DMSetApplicationContext(mshdm, &steel);  /* User Context */
>
>     KSPCreate(PETSC_COMM_WORLD, &ksp);
>     KSPSetDM(ksp, mshdm);
>     KSPSetComputeRHS(ksp, RHSFun, &steel);
>     KSPSetComputeOperators(ksp, JACFun, &steel);
>     KSPSetFromOptions(ksp);
>
> I have no problems with the ComputeRHS function that I wrote, it seems 
> to work fine. The ComputeOperators function on the other hand is 
> giving me an "Argument out of range error" when MatSetValues() is 
> called the second time (for e=1, discovered by debugging). Here's the 
> function definition:
>
>     PetscErrorCode JACFun(KSP ksp, Mat J, Mat jac, void* ctx)
>     {
>       Mater* pblm = (Mater*)ctx;
>       PetscInt N, ix, n;
>       PetscScalar Le;
>       PetscErrorCode ierr;
>       MatNullSpace nullspc;
>
>       PetscScalar idx[2], idy[2], elstiff[2][2];// idx, idy are ids
>     for MatSetValues()
>       DM mshdm;
>
>       ierr = KSPGetDM(ksp, &mshdm);
>       ierr = DMDAGetInfo(mshdm, NULL, &N, NULL, NULL, &n, NULL, NULL,
>     NULL, NULL, NULL, NULL, NULL, NULL);
>       DMDAGetCorners(mshdm, &ix, NULL, NULL, &n, NULL, NULL);
>
>       // Construct Element Stiffness
>       Le = pblm->L/N;
>       elstiff[0][0] = pblm->A*pblm->E/Le;
>       elstiff[0][1] = -pblm->A*pblm->E/Le;
>       elstiff[1][0] = -pblm->A*pblm->E/Le;
>       elstiff[1][1] = pblm->A*pblm->E/Le;
>
>     PetscInt is, ie;
>       MatGetOwnershipRange(jac, &is, &ie);
>       for (int e=is; e<((ie<N-1)?ie:(N-1)); ++e) {// Only looping over
>     local elements
>         idx[0] = e; idx[1] = e+1;
>         idy[0] = e; idy[1] = e+1;
>
>         // This is the version I want to get working, but throws an
>     argument out of range error
>         MatSetValuesBlocked(jac, 2, (const PetscInt*)idx, 2, (const
>     PetscInt*)idy,
>                     (const PetscScalar*)elstiff, ADD_VALUES);
>
>         // This version seemed to not construct the correct matrix
>     (overlapped elements were added)
>         /* MatSetValue(jac, e, e, elstiff[0][0], ADD_VALUES); */
>         /* MatSetValue(jac, e, e+1, elstiff[0][1], ADD_VALUES); */
>         /* MatSetValue(jac, e+1, e, elstiff[1][0], ADD_VALUES); */
>         /* MatSetValue(jac, e+1, e+1, elstiff[1][1], ADD_VALUES); */
>       }
>       MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
>       ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
>
>       // I'm doing a "free-free" simulation, so constraining out the
>     nullspace
>     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nullspc);
>       MatSetNullSpace(J, nullspc);
>       MatNullSpaceDestroy(&nullspc);
>       PetscFunctionReturn(0);
>     }
>
>
> Thank you,
> Nidish
-- 
Nidish
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200809/cf16d9bb/attachment-0001.html>