Good Morning,
I'm trying to verify that the CAD -> PETSc/DMPlex methods I've developed
can be used for FEM analyses using PETSc. Attached is my current attempt
where I import a CAD STEP file to create a volumetric tetrahedral
discretization (DMPlex), designate boundary condition points using
DMLabels, and solve the Laplace problem (heat) with Dirichlet conditions on
each end. At command line I indicate the STEP file with the -filename
option and the dual space degree with -petscspace_degree 2. The run ends
with either a SEGV Fault or a General MPI Communication Error.
Could you please look over the attached file to tell me if what I'm doing
to set up the FEM problem is wrong?
Thank you in advance for your time and help.
-Brandon
TYPICAL ERROR MESSAGE
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: General MPI error
[0]PETSC ERROR: MPI error 605109765 Invalid communicator, error stack:
PMPI_Comm_get_attr(344): MPI_Comm_get_attr(comm=0x0,
comm_keyval=-1539309568, attribute_val=0x7ffe75a58848, flag=0x7ffe75a58844)
failed
MPII_Comm_get_attr(257): MPIR_Comm_get_attr(comm=0x0,
comm_keyval=-1539309568, attribute_val=0x7ffe75a58848, flag=0x7ffe75a58844)
failed
MPII_Comm_get_attr(53).: Invalid communicator
[0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could
be the program crashed before they were used or a spelling mistake, etc!
[0]PETSC ERROR: Option left: name:-dm_plex_refine_without_snap_to_geom
value: 0 source: command line
[0]PETSC ERROR: Option left: name:-dm_refine value: 1 source: command line
[0]PETSC ERROR: Option left: name:-snes_monitor (no value) source:
command line
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.18.5-1817-gd2497b8de4c
GIT Date: 2023-05-22 18:44:03 +0000
[0]PETSC ERROR: ./thermal on a named XPS. by bdenton Wed Jun 7 11:03:43
2023
[0]PETSC ERROR: Configure options --with-make-np=16
--prefix=/mnt/c/Users/Brandon/software/libs/petsc/3.19.1-gitlab/gcc/11.2.0/mpich/3.4.2/openblas/0.3.17/opt
--with-debugging=false --COPTFLAGS="-O3 -mavx" --CXXOPTFLAGS="-O3 -mavx"
--FOPTFLAGS=-O3 --with-shared-libraries=1
--with-mpi-dir=/mnt/c/Users/Brandon/software/libs/mpich/3.4.2/gcc/11.2.0
--with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1
--with-parmetis=true --download-parmetis=1 --with-superlu=true
--download-superlu=1 --with-superludir=true --download-superlu_dist=1
--with-blacs=true --download-blacs=1 --with-scalapack=true
--download-scalapack=1 --with-hypre=true --download-hypre=1
--with-hdf5-dir=/mnt/c/Users/Brandon/software/libs/hdf5/1.12.1/gcc/11.2.0
--with-valgrind-dir=/mnt/c/Users/Brandon/software/apps/valgrind/3.14.0
--with-blas-lib="[/mnt/c/Users/Brandon/software/libs/openblas/0.3.17/gcc/11.2.0/lib/libopenblas.so]"
--with-lapack-lib="[/mnt/c/Users/Brandon/software/libs/openblas/0.3.17/gcc/11.2.0/lib/libopenblas.so]"
--LDFLAGS= --with-tetgen=true --download-tetgen=1 --download-ctetgen=1
--download-opencascade=1 --download-egads
[0]PETSC ERROR: #1 PetscObjectName() at
/mnt/c/Users/Brandon/software/builddir/petsc-3.19.1-gitlab/src/sys/objects/pname.c:119
[0]PETSC ERROR: #2 PetscObjectGetName() at
/mnt/c/Users/Brandon/software/builddir/petsc-3.19.1-gitlab/src/sys/objects/pgname.c:27
[0]PETSC ERROR: #3 PetscDSAddBoundary() at
/mnt/c/Users/Brandon/software/builddir/petsc-3.19.1-gitlab/src/dm/dt/interface/dtds.c:3404
[0]PETSC ERROR: #4 DMAddBoundary() at
/mnt/c/Users/Brandon/software/builddir/petsc-3.19.1-gitlab/src/dm/interface/dm.c:7828
[0]PETSC ERROR: #5 main() at
/mnt/c/Users/Brandon/Documents/School/Dissertation/Software/EGADS-dev/thermal_v319/thermal_nozzle.c:173
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -dm_plex_geom_print_model 1 (source: command line)
[0]PETSC ERROR: -dm_plex_geom_shape_opt 0 (source: command line)
[0]PETSC ERROR: -dm_plex_refine_without_snap_to_geom 0 (source: command
line)
[0]PETSC ERROR: -dm_refine 1 (source: command line)
[0]PETSC ERROR: -filename ./examples/Nozzle_example.stp (source: command
line)
[0]PETSC ERROR: -petscspace_degree 2 (source: command line)
[0]PETSC ERROR: -snes_monitor (source: command line)
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-ma...@mcs.anl.gov----------
application called MPI_Abort(MPI_COMM_SELF, 98) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=98
:
system msg for write_line failure : Bad file descriptor
static const char help[] = "Test of FEA DMPlex w/CAD functionality";
#include <petsc.h>
/* User-defined Work Context */
typedef struct {
PetscInt dummy;
char filename[PETSC_MAX_PATH_LEN]; // context containing CAD filename
from command line
} AppCtx;
void f0_function(PetscInt dim, PetscInt Nf, PetscInt NfAux,
const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar
u[], const PetscScalar u_t[], const PetscScalar u_x[],
const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar
a[], const PetscScalar a_t[], const PetscScalar a_x[],
PetscReal t, const PetscReal x[], PetscScalar f0[])
{
f0[0] = 0.0;
}
void f1_function(PetscInt dim, PetscInt Nf, PetscInt NfAux,
const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar
u[], const PetscScalar u_t[], const PetscScalar u_x[],
const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar
a[], const PetscScalar a_t[], const PetscScalar a_x[],
PetscReal t, const PetscReal x[], PetscScalar f1[])
{
for(PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d];
}
void g3_uu_function(PetscInt dim, PetscInt Nf, PetscInt NfAux,
const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar
u[], const PetscScalar u_t[], const PetscScalar u_x[],
const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar
a[], const PetscScalar a_t[], const PetscScalar a_x[],
PetscReal t, const PetscReal x[], PetscScalar g3[])
{
for (PetscInt d=0; d < dim; ++d) g3[d*dim + d] = 1.0;
}
//void bc_inlet(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc,
PetscScalar bcval[])
static PetscErrorCode bc_inlet(PetscInt dim, PetscReal time, const PetscReal
x[], PetscInt Nc, PetscScalar *u, void *ctx)
{
u[0] = 1400.0;
return PETSC_SUCCESS;
//bcval[0] = 1400.0;
}
//void bc_outlet(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt
Nc, PetscScalar bcval[])
static PetscErrorCode bc_outlet(PetscInt dim, PetscReal time, const PetscReal
x[], PetscInt Nc, PetscScalar *u, void *ctx)
{
u[0] = 100.0;
return PETSC_SUCCESS;
//bcval[0] = 100.0;
}
/* Procees Options - This should be removed once creation bug is fixed by Prof.
Knepley */
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
PetscFunctionBeginUser;
options->filename[0] = '\0';
PetscOptionsBegin(comm, "", "FEA DMPlex w/CAD Options", "DMPlex w/CAD");
PetscOptionsString("-filename", "The CAD/Geometry file", "ex18.c",
options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
PetscOptionsEnd();
PetscFunctionReturn(0);
}
/* Main Function */
int main(int argc, char **argv)
{
// Define PETSc Variables
DM dm, dmSurface; // Unstructured Grid
SNES snes; // Nonlinear Solver
Vec temp; // Solutions
AppCtx ctx; // User-defined Work Context
PetscInt dim; // DM Geometric Dimension
PetscBool simplex; // Is DMPlex Simplex type?
PetscFE fe; // PETSc Finite Element Object
PetscDS ds; // PETSc Discretization System Object
PetscMPIInt rank; // PETSc MPI Processor ID
PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &argv, NULL, help)); // Initialize PETSc
PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); // Process Options ::
Here we get the filename from the command line options
PetscCall(DMCreate(PETSC_COMM_WORLD, &dmSurface));
PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
PetscCall(DMSetType(dmSurface, DMPLEX));
PetscCall(DMSetType(dm, DMPLEX));
// Create DMPlex from CAD file
PetscCall(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); // Get Rank of
Current Processor
if (!rank) {
PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ctx.filename, "EGADS",
PETSC_TRUE, &dmSurface)); // Return dm created from CAD file
}
// Setup DMLabels for Boundary Conditions on DMPlex Vertices
PetscInt nStart, nEnd, eStart, eEnd, faceID;
DMLabel inletLabel, outletLabel;
const PetscInt idA = 1, idB = 2;
PetscCall(DMCreateLabel(dmSurface, "inlet"));
PetscCall(DMCreateLabel(dmSurface, "outlet"));
PetscCall(DMGetLabel(dmSurface, "inlet", &inletLabel));
PetscCall(DMGetLabel(dmSurface, "outlet", &outletLabel));
PetscCall(DMPlexGetDepthStratum(dmSurface, 0, &nStart, &nEnd));
PetscCall(DMPlexGetDepthStratum(dmSurface, 1, &eStart, &eEnd));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [nStart = %d || nEnd = %d] \n",
nStart, nEnd));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [eStart = %d || eEnd = %d] \n",
eStart, eEnd));
for (PetscInt ii = nStart; ii < nEnd; ++ii) {
PetscCall(DMGetLabelValue(dmSurface, "EGADS Face ID", ii, &faceID));
if (faceID == 14) {PetscCall(DMSetLabelValue(dmSurface, "inlet", ii,
faceID));}
if (faceID == 7) {PetscCall(DMSetLabelValue(dmSurface, "outlet", ii,
faceID));}
}
for (PetscInt ii = eStart; ii < eEnd; ++ii) {
PetscCall(DMGetLabelValue(dmSurface, "EGADS Face ID", ii, &faceID));
if (faceID == 14) {PetscCall(DMSetLabelValue(dmSurface, "inlet", ii,
faceID));}
if (faceID == 7) {PetscCall(DMSetLabelValue(dmSurface, "outlet", ii,
faceID));}
}
// Generate Volumetric Mesh
PetscCall(DMPlexGenerate(dmSurface, "tetgen", PETSC_TRUE, &dm)); //
Generate Volumetric Mesh
PetscCall(DMDestroy(&dmSurface)); // Destroy DM dmSurface no longer
needed
PetscCall(DMSetApplicationContext(dm, &ctx)); // Link context to dm
PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); // Write DM to file
(hdf5)
PetscCall(DMView(dm, PETSC_VIEWER_STDOUT_SELF));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [CREATED MESH] \n"));
// Setup SNES and link to DM
PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); // Create SNES object
PetscCall(SNESSetDM(snes, dm)); // Link SNES object
to DM
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SETUP SNES] \n"));
// Setup Discretization
PetscCall(DMGetDimension(dm, &dim)); // Get Geometric Dimension of the
DM (2D or 3D)
PetscCall(DMPlexIsSimplex(dm, &simplex)); // Check if DMPlex is Simplex Type
PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1,
simplex, NULL, -1, &fe)); // Define and Create PetscFE Object
PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); // Add Field to
Discretization object. In this case, the temperature field.
PetscCall(DMCreateDS(dm)); // Create The Discretization
System (DS) based on the field(s) added to the DM
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SETUP DISCRETIZATION] \n"));
// Setup Residuals & (Optional) Jacobian
PetscInt testFieldID = 0; // This is an assumption. How do
we find out which field ID goes with each added field?
PetscCall(DMGetDS(dm, &ds)); // Get DS associated with DM
PetscCall(PetscDSSetResidual(ds, testFieldID, (void (*))f0_function, (void
(*))f1_function)); // Set Residual Function
PetscCall(PetscDSSetJacobian(ds, testFieldID, testFieldID, NULL, NULL, NULL,
(void (*))g3_uu_function)); // Set Jacobian Function
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SETUP RESIDUAL & JACOBIAN] \n"));
// Apply inlet face conditions - Dirichlet
PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "inTemp", inletLabel, 1, &idA,
0, 0, NULL, (void (*)(void))bc_inlet, NULL, &ctx, NULL));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SETUP BOUNDARY CONDITION inTemp]
\n"));
PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "outTemp", outletLabel, 1, &idB,
0, 0, NULL, (void (*)(void))bc_outlet, NULL, &ctx, NULL));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SETUP BOUNDARY CONDITIONS] \n"));
// Create Global Vector, Initialize Values & Name
PetscCall(DMCreateGlobalVector(dm, &temp));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [AFTER GLOBAL VECTOR] \n"));
PetscCall(VecSet(temp, 100.0));
PetscCall(VecView(temp, PETSC_VIEWER_STDOUT_SELF));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [AFTER temp set to 100] \n"));
PetscCall(PetscObjectSetName((PetscObject)temp, "temperature"));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [AFTER SET OBJECT NAME] \n"));
PetscCall(DMPlexSetSNESLocalFEM(dm, &ctx, &ctx, &ctx));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [AFTER SNESLocalFEM] \n"));
//Set SNES from Options
PetscCall(SNESSetFromOptions(snes));
PetscCall(SNESView(snes, PETSC_VIEWER_STDOUT_SELF));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [AFTER SNESSetFromOptions] \n"));
// Solve System of Equations
PetscCall(SNESSolve(snes, NULL, temp));
PetscCall(PetscPrintf(PETSC_COMM_SELF, " [SOLVED PROBLEM] \n"));
//Get Solution and View
PetscCall(SNESGetSolution(snes, &temp));
PetscCall(VecViewFromOptions(temp, NULL, "-sol_view"));
// Cleanup
PetscCall(VecDestroy(&temp));
PetscCall(SNESDestroy(&snes));
PetscCall(DMDestroy(&dm));
PetscCall(PetscFinalize());
return 0;
}
