[petsc-users] Question about how to use DS to do the discretization

[Matthew Knepley]

On Sat, Mar 23, 2024 at 8:03 AM Gong Yujie <yc17470 at connect.um.edu.mo>
wrote:

> Dear PETSc group, I'm reading the DS part for the discretization start
> from SNES ex17. c which is a demo for solving linear elasticity problem. I
> have two questions for the details. The first question is for the residual
> function. Is the residual
> ZjQcmQRYFpfptBannerStart
> This Message Is From an External Sender
> This message came from outside your organization.
>
> ZjQcmQRYFpfptBannerEnd
> Dear PETSc group,
>
> I'm reading the DS part for the discretization start from SNES ex17.c
> which is a demo for solving linear elasticity problem. I have two questions
> for the details.
>
> The first question is for the residual function. Is the residual
> calculated as this? The dot product is a little weird because of the
> dimension of the result.
> Here \sigma is the stress tensor, \phi_i is the test function for the i-th
> function (Linear elasticity in 3D contains three equations).
>

The stress term in the momentum equation is

  (-div sigma) . psi = d_i sigma_{ij} psi_j

which is then integrated by parts

  sigma_{ij} d_i psi_j

This is linear isotropic elasticity, so

  sigma = \mu (d_i u_j + d_j u_i) + \lambda \delta_{ij} sigma_{kk}

In PETSc, we phrase the term in the weak form as

  f^1_{ij} d_i psi_j

so f1[i * dim + j] below is sigma_{ij}, from line 298 of ex17.c

  for (PetscInt c = 0; c < Nc; ++c) {
    for (PetscInt d = 0; d < dim; ++d) {
      f1[c * dim + d] += mu * (u_x[c * dim + d] + u_x[d * dim + c]);
      f1[c * dim + c] += lambda * u_x[d * dim + d];
    }
  }


> The second question is how to derive the Jacobian of the system (line 330
> in ex17.c). As shown in the PetscDSSetJacobian, we need to provide function
> g3() which I think is a 4th-order tensor with size 3*3*3*3 in this linear
> elasticity case. I'm not sure how to get it. Are there any references on
> how to get this Jacobian?
>

The Jacobian indices are shown here:
https://urldefense.us/v3/__https://petsc.org/main/manualpages/FE/PetscFEIntegrateJacobian/__;!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc53QKQec$ 
where the g3 term is

  \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u)
\nabla\phi^{gc}_g(q)

To get the Jacobian, we use u = \sum_i u_i psi_i, where psi_i is a vector,
and then differentiate the expression with respect to the coefficient u_i.
Since the operator is already linear,this is just matching indices

  for (PetscInt c = 0; c < Nc; ++c) {
    for (PetscInt d = 0; d < dim; ++d) {
      g3[((c * Nc + c) * dim + d) * dim + d] += mu;
      g3[((c * Nc + d) * dim + d) * dim + c] += mu;
      g3[((c * Nc + d) * dim + c) * dim + d] += lambda;
    }
  }

Take the first mu term, mu (d_c u_d ) (d_c psi_d). We know that fc == gc
and df == dg, so we get

      g3[((c * Nc + c) * dim + d) * dim + d] += mu;

For the second term, we transpose grad u, so fc == dg and gc == df,

      g3[((c * Nc + d) * dim + d) * dim + c] += mu;

Finally, for the lambda term, fc == df and gc == dg because we are matching
terms in the sum

      g3[((c * Nc + d) * dim + c) * dim + d] += lambda;

  Thanks,

     Matt


> I've checked about the comment before this Jacobian function (line 330 in
> ex17.c) but don't know how to get this.
>
> Thanks in advance!
>
> Best Regards,
> Yujie
>


-- 
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

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc4DkrV1c$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGcwApR5ek$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240323/9d7d5b3a/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 34587 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240323/9d7d5b3a/attachment-0001.png>