Hi Matthew
It seems like the problem is not fully fixed. I have changed the code to now
run on with both 2,3 and 4 cells. When I run the code using NP = 1..3 I get
different result both for NP=1 to NP=2/3 when cell count is larger than 2.
Kind regards,
Morten
____
mpiexec -np 1 ./ex18k
cells 2
Loc size: 36
Trace of matrix: 132.000000
cells 3
Loc size: 48
Trace of matrix: 192.000000
cells 4
Loc size: 60
Trace of matrix: 258.000000
mpiexec -np 2 ./ex18k
cells 2
Loc size: 24
Loc size: 24
Trace of matrix: 132.000000
cells 3
Loc size: 36
Loc size: 24
Trace of matrix: 198.000000
cells 4
Loc size: 36
Loc size: 36
Trace of matrix: 264.000000
mpiexec -np 3 ./ex18k
cells 2
Loc size: 24
Loc size: 24
Loc size: 0
Trace of matrix: 132.000000
cells 3
Loc size: 24
Loc size: 24
Loc size: 24
Trace of matrix: 198.000000
cells 4
Loc size: 36
Loc size: 24
Loc size: 24
Trace of matrix: 264.000000
________________________________
From: [email protected] [[email protected]] on
behalf of Morten Nobel-Jørgensen [[email protected]]
Sent: Sunday, September 25, 2016 11:15 AM
To: Matthew Knepley
Cc: PETSc [[email protected]]
Subject: Re: [petsc-users] DMPlex problem
Hi Matthew
Thank you for the bug-fix :) I can confirm that it works :)
And thanks for your hard work on PETSc - your work is very much appreciated!
Kind regards,
Morten
________________________________
From: Matthew Knepley [[email protected]]
Sent: Friday, September 23, 2016 2:46 PM
To: Morten Nobel-Jørgensen
Cc: PETSc [[email protected]]
Subject: Re: [petsc-users] DMPlex problem
On Fri, Sep 23, 2016 at 7:45 AM, Matthew Knepley
<[email protected]<mailto:[email protected]>> wrote:
On Fri, Sep 23, 2016 at 3:48 AM, Morten Nobel-Jørgensen
<[email protected]<mailto:[email protected]>> wrote:
Dear PETSc developers
Any update on this issue regarding DMPlex? Or is there any obvious workaround
that we are unaware of?
I have fixed this bug. It did not come up in nightly tests because we are not
using MatSetValuesLocal(). Instead we
use MatSetValuesClosure() which translates differently.
Here is the branch
https://bitbucket.org/petsc/petsc/branch/knepley/fix-dm-ltog-bs
and I have merged it to next. It will go to master in a day or two.
Also, here is the cleaned up source with no memory leaks.
Matt
Also should we additionally register the issue on Bitbucket or is reporting the
issue on the mailing list enough?
Normally we are faster, but the start of the semester was hard this year.
Thanks,
Matt
Kind regards,
Morten
________________________________
From: Matthew Knepley [[email protected]<mailto:[email protected]>]
Sent: Friday, September 09, 2016 12:21 PM
To: Morten Nobel-Jørgensen
Cc: PETSc [[email protected]<mailto:[email protected]>]
Subject: Re: [petsc-users] DMPlex problem
On Fri, Sep 9, 2016 at 4:04 AM, Morten Nobel-Jørgensen
<[email protected]<mailto:[email protected]>> wrote:
Dear PETSc developers and users,
Last week we posted a question regarding an error with DMPlex and multiple dofs
and have not gotten any feedback yet. This is uncharted waters for us, since we
have gotten used to an extremely fast feedback from the PETSc crew. So - with
the chance of sounding impatient and ungrateful - we would like to hear if
anybody has any ideas that could point us in the right direction?
This is my fault. You have not gotten a response because everyone else was
waiting for me, and I have been
slow because I just moved houses at the same time as term started here. Sorry
about that.
The example ran for me and I saw your problem. The local-tp-global map is
missing for some reason.
I am tracking it down now. It should be made by DMCreateMatrix(), so this is
mysterious. I hope to have
this fixed by early next week.
Thanks,
Matt
We have created a small example problem that demonstrates the error in the
matrix assembly.
Thanks,
Morten
--
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
--
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
--
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
static char help[] = "";
#include <petscdmplex.h>
const int dof = 3;
#undef __FUNCT__
#define __FUNCT__ "SetupDOFs"
/* Assign dofs (3 for each node) */
PetscErrorCode SetupDOFs(DM dm) {
PetscSection s;
PetscInt pStart, pEnd, cStart, cEnd, c, fStart, fEnd, f, vStart, vEnd, v, eStart, eEnd, e;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s); CHKERRQ(ierr);
ierr = PetscSectionSetNumFields(s, 1); CHKERRQ(ierr);
ierr = DMPlexGetChart(dm, &pStart, &pEnd); CHKERRQ(ierr);
PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); CHKERRQ(ierr);
ierr = PetscSectionSetChart(s, pStart, pEnd); CHKERRQ(ierr);
for (v = vStart; v < vEnd; ++v) {
ierr = PetscSectionSetDof(s, v, dof); CHKERRQ(ierr);
ierr = PetscSectionSetFieldDof(s, v, 0, dof); CHKERRQ(ierr);
}
ierr = PetscSectionSetUp(s); CHKERRQ(ierr);
ierr = DMSetDefaultSection(dm, s); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "CreateGlobalStiffnessMatrix"
PetscInt CreateGlobalStiffnessMatrix(DM dm, Mat *K) {
PetscSection localSection;
PetscInt cStart, cEnd, nStart, nEnd, e;
PetscInt edof[dof * 8]; // edof vector
PetscScalar ke[(dof * 8) * (dof * 8)]; // element matrix
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
ierr = DMPlexGetDepthStratum(dm, 0, &nStart, &nEnd); CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, K); CHKERRQ(ierr);
ierr = MatViewFromOptions(*K, NULL, "-mat_view"); CHKERRQ(ierr);
ierr = DMGetDefaultSection(dm, &localSection); CHKERRQ(ierr);
for (e = 0; e < (8 * dof) * (8 * dof); e++) {
ke[e] = 1.0; // placeholder value - replace with local stiffness matrix
}
for (e = cStart; e < cEnd; e++) {
PetscBool useCone = PETSC_TRUE;
PetscInt clojureCount;
PetscInt *clojureIndex = NULL;
PetscInt edofIndex = 0, v;
ierr = DMPlexGetTransitiveClosure(dm, e, useCone, &clojureCount, &clojureIndex); CHKERRQ(ierr);
for (v = 0; v < clojureCount; ++v) {
PetscInt cPoint = clojureIndex[v * 2];
PetscInt cOrientation = clojureIndex[v * 2 + 1];
PetscBool isVertex = cPoint >= nStart && cPoint < nEnd ? PETSC_TRUE : PETSC_FALSE; // is vertex if it exist in the vertex range
if (isVertex) {
PetscInt offset, dofSize, i;
ierr = PetscSectionGetOffset(localSection, cPoint, &offset); CHKERRQ(ierr);
ierr = PetscSectionGetDof(localSection, cPoint, &dofSize); CHKERRQ(ierr);
for (i = 0; i < dofSize; ++i) {
edof[edofIndex] = offset + i;
edofIndex++;
}
}
}
ierr = DMPlexRestoreTransitiveClosure(dm, e, useCone, &clojureCount, &clojureIndex); CHKERRQ(ierr);
ierr = MatSetValuesLocal(*K, 8 * dof, edof, 8 * dof, edof, ke, ADD_VALUES); CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(*K, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(*K, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv) {
DM dm, dist;
Mat mat;
PetscInt dim = 3;
PetscInt overlap = 0;
PetscErrorCode ierr;
PetscViewer view;
ierr = PetscInitialize(&argc, &argv, NULL, help); CHKERRQ(ierr);
for (int c = 2;c <= 4;c++){
PetscPrintf(PETSC_COMM_WORLD, "cells %i\n", c);
PetscInt cells[] = {1, 1, c};
ierr = DMPlexCreateHexBoxMesh(PETSC_COMM_WORLD, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
&dm); CHKERRQ(ierr);
ierr = DMPlexSetAdjacencyUseCone(dm, PETSC_FALSE); CHKERRQ(ierr);
ierr = DMPlexSetAdjacencyUseClosure(dm, PETSC_TRUE); CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, NULL, &dist); CHKERRQ(ierr);
if (dist) {
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dist;
}
ierr = SetupDOFs(dm); CHKERRQ(ierr);
ierr = CreateGlobalStiffnessMatrix(dm, &mat); CHKERRQ(ierr);
Vec v;
ierr = DMCreateLocalVector(dm, &v); CHKERRQ(ierr);
PetscInt vsize;
ierr = VecGetSize(v, &vsize); CHKERRQ(ierr);
PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Loc size: %i\n", vsize);
PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
PetscScalar trace;
ierr = MatGetTrace(mat, &trace); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "Trace of matrix: %f\n", trace);
}
PetscFinalize();
return 0;
}
