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

[Nidish]

Update: It was a type error (idx and idy below were declared as 
PetscScalar* instead of PetscInt*), I'm sorry.

Guess I'm being punished for not paying attention to types.

Nidish

On 8/9/20 10:05 PM, Nidish wrote:
>
> 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
-- 
Nidish
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200809/f4d26564/attachment.html>