<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Update: It was a type error (idx and idy below were declared as
PetscScalar* instead of PetscInt*), I'm sorry. <br>
</p>
<p>Guess I'm being punished for not paying attention to types.</p>
<p>Nidish<br>
</p>
<div class="moz-cite-prefix">On 8/9/20 10:05 PM, Nidish wrote:<br>
</div>
<blockquote type="cite"
cite="mid:d94e0663-3157-5dab-1ba6-816382f98790@rice.edu">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<p>Update: I had made some mistakes in the assembly in the
previous version. They're fixed now (see below). <br>
</p>
<div class="moz-cite-prefix">On 8/9/20 9:13 PM, Nidish wrote:<br>
</div>
<blockquote type="cite"
cite="mid:4960b5fb-d8fc-ecbf-4083-618a8a8d6642@rice.edu">
<meta http-equiv="content-type" content="text/html;
charset=UTF-8">
<p>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.</p>
<p>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:</p>
<blockquote>
<pre><a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode" moz-do-not-send="true">PetscErrorCode</a> <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/MatSetValues.html#MatSetValues" moz-do-not-send="true">MatSetValues</a>(<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/Mat.html#Mat" moz-do-not-send="true">Mat</a> mat,<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> m,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> idxm[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> n,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> idxn[],const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscScalar.html#PetscScalar" moz-do-not-send="true">PetscScalar</a> v[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/InsertMode.html#InsertMode" moz-do-not-send="true">InsertMode</a> addv)
</pre>
</blockquote>
<div class="moz-signature">Here's the DM setup for my code: <tt><br>
</tt></div>
<blockquote>
<div class="moz-signature"><tt>DM mshdm;</tt><tt><br>
</tt><tt>Mater steel = {1.0, 2e11, 7800, 3.0, 400.0};</tt><tt>
// A structure with {Area, Modulus, Density, Length,
Force} (runtime vars)<br>
</tt><tt>KSP ksp;</tt><tt><br>
</tt><tt><br>
</tt><tt>DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N,
1, 1, NULL, &mshdm);</tt><tt> // N is a runtime
variable<br>
</tt><tt>DMSetFromOptions(mshdm);</tt><tt><br>
</tt><tt>DMSetUp(mshdm);</tt><tt><br>
</tt><tt>DMSetApplicationContext(mshdm, &steel); /*
User Context */</tt><tt><br>
</tt><tt><br>
</tt><tt>KSPCreate(PETSC_COMM_WORLD, &ksp);</tt><tt><br>
</tt><tt>KSPSetDM(ksp, mshdm);</tt><tt><br>
</tt><tt>KSPSetComputeRHS(ksp, RHSFun, &steel);</tt><tt><br>
</tt><tt>KSPSetComputeOperators(ksp, JACFun, &steel);</tt><tt><br>
</tt><tt>KSPSetFromOptions(ksp);</tt><br>
</div>
</blockquote>
<div class="moz-signature">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:</div>
<blockquote>
<div class="moz-signature"><tt>PetscErrorCode JACFun(KSP ksp,
Mat J, Mat jac, void* ctx)</tt><tt><br>
</tt><tt>{</tt><tt><br>
</tt><tt> Mater* pblm = (Mater*)ctx;</tt><tt><br>
</tt><tt> PetscInt N, ix, n;</tt><tt><br>
</tt><tt> PetscScalar Le;</tt><tt><br>
</tt><tt> PetscErrorCode ierr;</tt><tt><br>
</tt><tt> MatNullSpace nullspc;</tt><tt><br>
</tt><tt><br>
</tt><tt> PetscScalar idx[2], idy[2], elstiff[2][2];</tt><tt>
// idx, idy are ids for MatSetValues()<br>
</tt><tt> DM mshdm;</tt><tt><br>
</tt><tt><br>
</tt><tt> ierr = KSPGetDM(ksp, &mshdm);</tt><tt><br>
</tt><tt> ierr = DMDAGetInfo(mshdm, NULL, &N, NULL,
NULL, &n, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
NULL);</tt><tt><br>
</tt><tt> DMDAGetCorners(mshdm, &ix, NULL, NULL,
&n, NULL, NULL);</tt><tt><br>
</tt><tt><br>
</tt><tt> // Construct Element Stiffness</tt></div>
<div class="moz-signature"><tt> Le = pblm->L/N;</tt><tt><br>
</tt><tt> elstiff[0][0] = pblm->A*pblm->E/Le;</tt><tt><br>
</tt><tt> elstiff[0][1] = -pblm->A*pblm->E/Le;</tt><tt><br>
</tt><tt> elstiff[1][0] = -pblm->A*pblm->E/Le; </tt><tt><br>
</tt><tt> elstiff[1][1] = pblm->A*pblm->E/Le;<br>
</tt><br>
<tt>PetscInt is, ie;</tt><tt><br>
</tt><tt> MatGetOwnershipRange(jac, &is, &ie);</tt><tt><br>
</tt><tt> for (int e=is; e<((ie<N-1)?ie:(N-1)); ++e)
{</tt><tt> // Only looping over local elements<br>
</tt><tt> idx[0] = e; idx[1] = e+1;</tt><tt><br>
</tt><tt> idy[0] = e; idy[1] = e+1; </tt><tt><br>
</tt><tt><br>
</tt><tt> // This is the version I want to get working,
but throws an argument out of range error<br>
</tt></div>
<div class="moz-signature"><tt> MatSetValuesBlocked(jac, 2,
(const PetscInt*)idx, 2, (const PetscInt*)idy,</tt><tt><br>
</tt><tt> (const PetscScalar*)elstiff,
ADD_VALUES);</tt></div>
<div class="moz-signature"><tt><br>
</tt><tt> // This version seemed to not construct the
correct matrix (overlapped elements were added)<br>
</tt><tt> /* MatSetValue(jac, e, e, elstiff[0][0],
ADD_VALUES); */</tt><tt><br>
</tt><tt> /* MatSetValue(jac, e, e+1, elstiff[0][1],
ADD_VALUES); */</tt><tt><br>
</tt><tt> /* MatSetValue(jac, e+1, e, elstiff[1][0],
ADD_VALUES); */</tt><tt><br>
</tt><tt> /* MatSetValue(jac, e+1, e+1, elstiff[1][1],
ADD_VALUES); */</tt><tt><br>
</tt><tt> }</tt><tt><br>
</tt><tt> MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);</tt><tt><br>
</tt><tt> ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);</tt><tt><br>
</tt><tt><br>
</tt><tt> // I'm doing a "free-free" simulation, so
constraining out the nullspace</tt></div>
<div class="moz-signature"><tt>
MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL,
&nullspc);</tt><tt><br>
</tt><tt> MatSetNullSpace(J, nullspc);</tt><tt><br>
</tt><tt> MatNullSpaceDestroy(&nullspc);</tt><tt><br>
</tt><tt> PetscFunctionReturn(0);</tt><tt><br>
</tt><tt>}</tt><br>
</div>
</blockquote>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Thank you,<br>
Nidish</div>
</blockquote>
<div class="moz-signature">-- <br>
Nidish</div>
</blockquote>
<div class="moz-signature">-- <br>
Nidish</div>
</body>
</html>