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
