Re: [petsc-users] Fortran interfaces: Google Summer of Code project?
Looks good to me. Thanks for taking the initiative. Martin Diehl writes: > Dear All, > > pls. find attached the proposal for > https: //urldefense. us/v3/__https: //github. com/fortran-lang/webpage/wiki/GSoC-2024-Project-ideas__;!!G_uCfscf7eWS!blt-tme3gpPtGXW8o_NasIEtvlUql08d2_Q9Q83PvjPwLqIITWec4gsejnmPLpUCgMU8xuIEZgTlg65QIRo$. ZjQcmQRYFpfptBannerStart This Message Is From an External Sender This message came from outside your organization. ZjQcmQRYFpfptBannerEnd Looks good to me. Thanks for taking the initiative. Martin Diehl writes: > Dear All, > > pls. find attached the proposal for > https://urldefense.us/v3/__https://github.com/fortran-lang/webpage/wiki/GSoC-2024-Project-ideas__;!!G_uCfscf7eWS!blt-tme3gpPtGXW8o_NasIEtvlUql08d2_Q9Q83PvjPwLqIITWec4gsejnmPLpUCgMU8xuIEZgTlg65QIRo$. > > I tried to keep it general such that we keep all options open. > > Let me know if you want to change/add/remove anything and/or if you > want to be listed as mentor. > > Since I've mixed up the deadline, the most urgent task is the finding > of suitable candidates. Once it's online, I'll post on linkedin but > ideally we can motivate someone who is already known. > > best regards, > Martin > > On Thu, 2024-03-21 at 23:13 -0600, Jed Brown wrote: >> Barry Smith writes: >> >> > > We already have the generated ftn-auto-interfaces/*.h90. The >> > > INTERFACE keyword could be replaced with CONTAINS (making these >> > > definitions instead of just interfaces), and then the bodies >> > > could use iso_c_binding to call the C functions. That would >> > > reduce repetition and be the standards-compliant way to do this. >> > >> > Sure, the interface and the stub go in the same file instead of >> > two files. This is slightly nicer but not significantly simpler, >> > and alone, it is not reason enough to write an entire new stub >> > generator. >> >> I agree, but if one *is* writing a new stub generator for good >> reasons (like better automation/completeness), there's a case for >> doing it this way unless users really need an environment in which >> that feature can't be used. >> >> > > What we do now with detecting the Fortran mangling scheme and >> > > calling conventions "works" but doesn't conform to any standard >> > > and there's nothing stopping Fortran implementations from >> > > creating yet another variant that we have to deal with manually. >> > >> > From practical experience, calling C/Fortran using non-standards >> > has only gotten easier over the last thirty-five years (fewer >> > variants on how char* is handled); it has not gotten more >> > complicated, so I submit the likelihood of "nothing stopping >> > Fortran implementations from creating yet another variant that we >> > have to deal with manually" is (though possible) rather unlikely. >> > As far as I am concerned, much of iso_c_binding stuff just solved a >> > problem that never really existed (except in some people's minds) >> > since calling C/Fortran has always been easy, modulo knowing a tiny >> > bit of information.. >> >> An examples for concreteness: >> >> https://urldefense.us/v3/__https://fortranwiki.org/fortran/show/Generating*C*Interfaces__;Kys!!G_uCfscf7eWS!blt-tme3gpPtGXW8o_NasIEtvlUql08d2_Q9Q83PvjPwLqIITWec4gsejnmPLpUCgMU8xuIEZgTlsKy54Ps$ >> >> And discussion: >> >> https://urldefense.us/v3/__https://fortran-lang.discourse.group/t/iso-c-binding-looking-for-practical-example-of-how-it-helps-with-mangling/3393/8__;!!G_uCfscf7eWS!blt-tme3gpPtGXW8o_NasIEtvlUql08d2_Q9Q83PvjPwLqIITWec4gsejnmPLpUCgMU8xuIEZgTlVTiuUTk$ >> >> With this approach, one could even use method syntax like >> ksp%SetOperators(J, Jpre), as in the nlopt-f project linked in the >> top of this question. I don't know if we want that (it would be a >> huge change for users, albeit it "easy"), but generating stubs in >> Fortran using iso_c_binding opens a variety of possibilities for more >> idiomatic bindings. > > -- > KU Leuven > Department of Computer Science > Department of Materials Engineering > Celestijnenlaan 200a > 3001 Leuven, Belgium > # Improved generation of Fortran interfaces for [PETSc](https://urldefense.us/v3/__https://petsc.org__;!!G_uCfscf7eWS!blt-tme3gpPtGXW8o_NasIEtvlUql08d2_Q9Q83PvjPwLqIITWec4gsejnmPLpUCgMU8xuIEZgTljowbJv0$) > > PETSc, the Portable, Extensible Toolkit for Scientific Computation, pronounced PET-see, is for the scalable (parallel) solution of scientific applications modeled by partial differential equations (PDEs). > It has bindings for C, Fortran, and Python (via petsc4py). PETSc also contains TAO, the Toolkit for Advanced Optimization, software library. It supports MPI, and GPUs through CUDA, HIP, Kokkos, or OpenCL, as well as hybrid MPI-GPU parallelism; it also supports the NEC-SX Tsubasa Vector Engine. > > Currently, only a part of the Fortran interfaces can be generated automatically using
Re: [petsc-users] Question about how to use DS to do the discretization
On Sat, Mar 23, 2024 at 8:03 AM Gong Yujie
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$
Re: [petsc-users] Question about how to use DS to do the discretization
Take a look at https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex11.c?ref_type=heads__;!!G_uCfscf7eWS!er4CI8GIe7OCWvCmRKQpZt6FOz1QYvbuZOdf2Fm7pvMGee3I9M5bhjNytv42F9C17NpBy0i6mTfgEmQfUR_QOqwC7gC6pYk$ and the discussion at the beginning (including the reference to the original paper) On Sat, Mar 23, 2024, 15:03 Gong Yujie 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 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? > > 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 >
[petsc-users] Question about how to use DS to do the discretization
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. [cid:b790e49a-5adb-4085-807c-546e7c2d941f] 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 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? 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
