[petsc-users] Drawing a line plot with two lines
Hi All,
I am trying to plot a loglog plot with two lines, one for my error and one with
a slope of 1. Plotting only the error or only the slope of 1 works fine, but I
like both lines in the same plot (that doesn't work till now). Does somebody
know how to solve this? Find a code snipet below:
PetscDraw draw;
PetscDrawLG lg;
PetscDrawAxis axis;
PetscReal xc, yc;
PetscDrawCreate(PETSC_COMM_SELF, NULL, "Log(Error) vs Log(dx)", PETSC_DECIDE,
PETSC_DECIDE, PETSC_DRAW_HALF_SIZE, PETSC_DRAW_HALF_SIZE, );
PetscDrawSetFromOptions(draw);
PetscDrawLGCreate(draw, 2, );
PetscDrawLGSetUseMarkers(lg, PETSC_TRUE);
PetscDrawLGGetAxis(lg, );
PetscDrawAxisSetLabels(axis, NULL, "Log(dx)", "Log(Error)");
for loop {
xc = PetscLog10Real(dx);
yc = PetscLog10Real(error);
PetscDrawLGAddPoint(lg, , ); // to plot the error
PetscDrawLGAddPoint(lg, , ); // to plot line with slope 1
PetscDrawLGDraw(lg);
}
PetscDrawSetPause(draw, -2);
PetscDrawLGDestroy();
PetscDrawDestroy();
Best regards,
Thijs Smit
Re: [petsc-users] Outputting cell data in stead of point data while writing .vtr file
Hi Matt, Oke, would a code change be difficult? I mean, feasible for me to do as a mechanical engineer? ;) Or can I use an other version of PETSc where the output to ASCII VTK is still available? Best, Thijs From: Matthew Knepley Sent: 11 March 2021 14:54 To: Smit Thijs Cc: petsc-users@mcs.anl.gov Subject: Re: [petsc-users] Outputting cell data in stead of point data while writing .vtr file On Thu, Mar 11, 2021 at 8:17 AM Smit Thijs mailto:thijs.s...@hest.ethz.ch>> wrote: Hi Matt, Actually I have two 3D DMDA’s, one for the nodal data, where the FEM is solved on. The other DMDA is a cell centered one for the volume data, like the density of a particular voxel. Ideally I would like to write both point data (displacement field) and cell data (density) to the vtr. Code for DMDA. DMBoundaryType bx = DM_BOUNDARY_NONE; DMBoundaryType by = DM_BOUNDARY_NONE; DMBoundaryType bz = DM_BOUNDARY_NONE; DMDAStencilType stype = DMDA_STENCIL_BOX; PetscInt stencilwidth = 1; // Create the nodal mesh ierr = DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stype, nx, ny, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, numnodaldof, stencilwidth, 0, 0, 0, &(da_nodes)); CHKERRQ(ierr); DMSetFromOptions(da_nodes); DMSetUp(da_nodes); ierr = DMDASetUniformCoordinates(da_nodes, xmin, xmax, ymin, ymax, zmin, zmax); CHKERRQ(ierr); ierr = DMDASetElementType(da_nodes, DMDA_ELEMENT_Q1); CHKERRQ(ierr); When I wrote the original version which output to ASCII VTK, we allowed switching between point data and cell data. It is fragile since you have to assure that the different grids match properly. When the viewer was rewritten to use the XML, all output was point data. It looks like it would take code changes to get this done. Thanks, Matt Best, Thijs From: Matthew Knepley mailto:knep...@gmail.com>> Sent: 11 March 2021 14:08 To: Smit Thijs mailto:thijs.s...@hest.ethz.ch>> Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov> Subject: Re: [petsc-users] Outputting cell data in stead of point data while writing .vtr file What kind of DM is it? Thanks, Matt On Thu, Mar 11, 2021 at 3:36 AM Smit Thijs mailto:thijs.s...@hest.ethz.ch>> wrote: Hi All, I am outputting several vectors to a .vtr file successfully for viewing in Paraview. At this moment the information is written to point data. How can I change this and make sure the data is written to cell data? The code I am currently using for outputting: PetscViewer viewer; ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD, “test.vtr”, FILE_MODE_WRITE, ); CHKERRQ(ierr); ierr = DMView(nd, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)xPhys,"xPhys"); ierr = VecView(xPhys, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)S,"SvonMises"); ierr = VecView(S, viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(); CHKERRQ(ierr); Best regards, Thijs Smit PhD Candidate ETH Zurich Institute for Biomechanics -- 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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/> -- 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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
Re: [petsc-users] Outputting cell data in stead of point data while writing .vtr file
Hi Matt, Actually I have two 3D DMDA’s, one for the nodal data, where the FEM is solved on. The other DMDA is a cell centered one for the volume data, like the density of a particular voxel. Ideally I would like to write both point data (displacement field) and cell data (density) to the vtr. Code for DMDA. DMBoundaryType bx = DM_BOUNDARY_NONE; DMBoundaryType by = DM_BOUNDARY_NONE; DMBoundaryType bz = DM_BOUNDARY_NONE; DMDAStencilType stype = DMDA_STENCIL_BOX; PetscInt stencilwidth = 1; // Create the nodal mesh ierr = DMDACreate3d(PETSC_COMM_WORLD, bx, by, bz, stype, nx, ny, nz, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, numnodaldof, stencilwidth, 0, 0, 0, &(da_nodes)); CHKERRQ(ierr); DMSetFromOptions(da_nodes); DMSetUp(da_nodes); ierr = DMDASetUniformCoordinates(da_nodes, xmin, xmax, ymin, ymax, zmin, zmax); CHKERRQ(ierr); ierr = DMDASetElementType(da_nodes, DMDA_ELEMENT_Q1); CHKERRQ(ierr); Best, Thijs From: Matthew Knepley Sent: 11 March 2021 14:08 To: Smit Thijs Cc: petsc-users@mcs.anl.gov Subject: Re: [petsc-users] Outputting cell data in stead of point data while writing .vtr file What kind of DM is it? Thanks, Matt On Thu, Mar 11, 2021 at 3:36 AM Smit Thijs mailto:thijs.s...@hest.ethz.ch>> wrote: Hi All, I am outputting several vectors to a .vtr file successfully for viewing in Paraview. At this moment the information is written to point data. How can I change this and make sure the data is written to cell data? The code I am currently using for outputting: PetscViewer viewer; ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD, “test.vtr”, FILE_MODE_WRITE, ); CHKERRQ(ierr); ierr = DMView(nd, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)xPhys,"xPhys"); ierr = VecView(xPhys, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)S,"SvonMises"); ierr = VecView(S, viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(); CHKERRQ(ierr); Best regards, Thijs Smit PhD Candidate ETH Zurich Institute for Biomechanics -- 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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
[petsc-users] Outputting cell data in stead of point data while writing .vtr file
Hi All, I am outputting several vectors to a .vtr file successfully for viewing in Paraview. At this moment the information is written to point data. How can I change this and make sure the data is written to cell data? The code I am currently using for outputting: PetscViewer viewer; ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD, "test.vtr", FILE_MODE_WRITE, ); CHKERRQ(ierr); ierr = DMView(nd, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)xPhys,"xPhys"); ierr = VecView(xPhys, viewer); CHKERRQ(ierr); PetscObjectSetName((PetscObject)S,"SvonMises"); ierr = VecView(S, viewer); CHKERRQ(ierr); ierr = PetscViewerDestroy(); CHKERRQ(ierr); Best regards, Thijs Smit PhD Candidate ETH Zurich Institute for Biomechanics
Re: [petsc-users] Loading external array data into PETSc
Hi Matt,
Thank you, I needed to reconfigure with hdf5.
The import works well now.
Best, Thijs
From: Matthew Knepley
Sent: 04 March 2021 14:53
To: Smit Thijs
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] Loading external array data into PETSc
On Thu, Mar 4, 2021 at 1:10 AM Smit Thijs
mailto:thijs.s...@hest.ethz.ch>> wrote:
Hi Matt, Roland and others,
@ Roland, I can create the Vec in Petsc, but want to load it with external
data. I can arrange that data in what ever way with Python. It is about two
vectors, each of length 4.6 x10^6 entries.
I have created a hdf5 file with Python and try to load it into PETSc which
gives me some errors. I added the small test script I am using. The error I am
getting is:
error: ‘PetscViewerHDF5Open’ was not declared in this scope; did you mean
‘PetscViewerVTKOpen’?
Am I using the wrong header
Did you configure PETSc with HDF5? If so, the configure.log will have an entry
for it at the bottom. If not, reconfigure with it:
${PETSC_DIR}/${PETSC_ARCH}/lib/petsc/conf/reconfigure-${PETSC_ARCH}.py
--download-hdf5
or you can use
--with-hdf5-dir=/path/to/hdf5
if it is already installed.
Thanks,
Matt
Best,
Thijs
From: Matthew Knepley mailto:knep...@gmail.com>>
Sent: 03 March 2021 23:02
To: Smit Thijs mailto:thijs.s...@hest.ethz.ch>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] Loading external array data into PETSc
On Wed, Mar 3, 2021 at 4:23 PM Smit Thijs
mailto:thijs.s...@hest.ethz.ch>> wrote:
Hi All,
I would like to readin a fairly large external vector into PETSc to be used
further. My idea was to write a hdf5 file using Python with the vector data.
Then read that hdf5 file into PETSc using PetscViewerHDF5Open and VecLoad. Is
this the advised way or is there a better alternative to achieve the same
(getting this external array data into PETSc)?
I think there are at least two nice ways to do this:
1) Use HDF5. Here you have to name your vector to match the HDF5 object
2) Use raw binary. This is a specific PETSc format, but it is fast and simple.
Thanks,
Matt
Best regards,
Thijs Smit
PhD Candidate
ETH Zurich
Institute for Biomechanics
--
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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
--
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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
Re: [petsc-users] Loading external array data into PETSc
Hi Matt, Roland and others,
@ Roland, I can create the Vec in Petsc, but want to load it with external
data. I can arrange that data in what ever way with Python. It is about two
vectors, each of length 4.6 x10^6 entries.
I have created a hdf5 file with Python and try to load it into PETSc which
gives me some errors. I added the small test script I am using. The error I am
getting is:
error: ‘PetscViewerHDF5Open’ was not declared in this scope; did you mean
‘PetscViewerVTKOpen’?
Am I using the wrong header?
Best,
Thijs
From: Matthew Knepley
Sent: 03 March 2021 23:02
To: Smit Thijs
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] Loading external array data into PETSc
On Wed, Mar 3, 2021 at 4:23 PM Smit Thijs
mailto:thijs.s...@hest.ethz.ch>> wrote:
Hi All,
I would like to readin a fairly large external vector into PETSc to be used
further. My idea was to write a hdf5 file using Python with the vector data.
Then read that hdf5 file into PETSc using PetscViewerHDF5Open and VecLoad. Is
this the advised way or is there a better alternative to achieve the same
(getting this external array data into PETSc)?
I think there are at least two nice ways to do this:
1) Use HDF5. Here you have to name your vector to match the HDF5 object
2) Use raw binary. This is a specific PETSc format, but it is fast and simple.
Thanks,
Matt
Best regards,
Thijs Smit
PhD Candidate
ETH Zurich
Institute for Biomechanics
--
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://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
#include
#include
#include
static char help[] = "Test HDF5\n\n";
int main(int argc, char* argv[]) {
// Error code for debugging
PetscErrorCode ierr;
// Initialize PETSc / MPI and pass input arguments to PETSc
PetscInitialize(, , PETSC_NULL, help);
Vec xPassive;
VecCreate(PETSC_COMM_WORLD, );
PetscObjectSetName((PetscObject) xPassive, "xPassive");
VecSetSizes(xPassive, PETSC_DECIDE, 4614720);
VecSetFromOptions(xPassive);
PetscViewer viewin;
PetscViewerHDF5Open(PETSC_COMM_WORLD, "xPassive.h5", FILE_MODE_READ,
);
PetscViewerDestroy();
VecDestroy();
PetscFinalize();
return 0;
}
[petsc-users] Loading external array data into PETSc
Hi All, I would like to readin a fairly large external vector into PETSc to be used further. My idea was to write a hdf5 file using Python with the vector data. Then read that hdf5 file into PETSc using PetscViewerHDF5Open and VecLoad. Is this the advised way or is there a better alternative to achieve the same (getting this external array data into PETSc)? Best regards, Thijs Smit PhD Candidate ETH Zurich Institute for Biomechanics
[petsc-users] error: invalid types ‘PetscScalar* {aka double*}[PetscScalar {aka double}]’ for array subscript
Hi All,
I am having the following error when I try to do a mapping with vectors and I
can’t figure out how to solve this or what is going wrong:
error: invalid types ‘PetscScalar* {aka double*}[PetscScalar {aka double}]’ for
array subscript
xpMMA[i] = xp[indicesMap[i]];
Herewith two code snippets:
// total number of elements on core
PetscInt nel;
VecGetLocalSize(xPhys, );
// create xPassive vector
ierr = VecDuplicate(xPhys, );
CHKERRQ(ierr);
// create mapping vector
ierr = VecDuplicate(xPhys, );
CHKERRQ(ierr);
// index set for xPassive and indicator
PetscScalar *xpPassive, *xpIndicator;
ierr = VecGetArray(xPassive, );
CHKERRQ(ierr);
ierr = VecGetArray(indicator, );
CHKERRQ(ierr);
// counters for total and active elements on this processor
PetscInt tcount = 0; // total number of elements
PetscInt acount = 0; // number of active elements
PetscInt scount = 0; // number of solid elements
PetscInt rcount = 0; // number of rigid element
// loop over all elements and update xPassive from wrapper data
// count number of active elements, acount
// set indicator vector
for (PetscInt el = 0; el < nel; el++) {
if (data.xPassive_w.size() > 1) {
xpPassive[el] = data.xPassive_w[el];
tcount++;
if (xpPassive[el] < 0) {
xpIndicator[acount] = el;
acount++;
}
} else {
xpPassive[el] = -1.0; // default, if no xPassive_w than all
elements are active = -1.0
}
}
// printing
//PetscPrintf(PETSC_COMM_WORLD, "tcount: %i\n", tcount);
//PetscPrintf(PETSC_COMM_WORLD, "acount: %i\n", acount);
// Allreduce, get number of active elements over all processes
// tmp number of var on proces
// acount total number of var sumed
PetscInt tmp = acount;
acount = 0.0;
MPI_Allreduce(, &(acount), 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD);
create xMMA vector
VecCreateMPI(PETSC_COMM_WORLD, tmp, acount, );
// Pointers to the vectors
PetscScalar *xp, *xpMMA, *indicesMap;
//PetscInt indicesMap;
ierr = VecGetArray(MMAVector, );
CHKERRQ(ierr);
ierr = VecGetArray(elementVector, );
CHKERRQ(ierr);
// Index set
PetscInt nLocalVar;
VecGetLocalSize(xMMA, );
// print number of var on pocessor
PetscPrintf(PETSC_COMM_WORLD, "Local var: %i\n", nLocalVar);
ierr = VecGetArray(indicator, );
CHKERRQ(ierr);
// Run through the indices
for (PetscInt i = 0; i < nLocalVar; i++) {
if (updateDirection > 0) {
//PetscPrintf(PETSC_COMM_WORLD, "i: %i, xp[%i] = %f\n", i,
indicesMap[i], xp[indicesMap[i]]);
xpMMA[i] = xp[indicesMap[i]];
} else if (updateDirection < 0) {
xp[indicesMap[i]] = xpMMA[i];
//PetscPrintf(PETSC_COMM_WORLD, "i: %i, xp[%i] = %f\n", i,
indicesMap[i], xp[indicesMap[i]]);
}
}
// Restore
ierr = VecRestoreArray(elementVector, );
CHKERRQ(ierr);
ierr = VecRestoreArray(MMAVector, );
CHKERRQ(ierr);
ierr = VecRestoreArray(indicator, );
CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "FINISHED UpdateVariables \n");
The error message says that the type with which I try to index is wrong, I
think. But VecGetArray only excepts scalars. Furthermore, the el variable is an
int, but is seams like to turn out to be a scalar. Does anybody see how to
proceed with this?
Best regards,
Thijs Smit
