[petsc-users] Modeling a Network with DMPlex?
Matthew Knepley
knepley at gmail.com
Mon Feb 17 08:39:58 CST 2014
On Mon, Feb 17, 2014 at 8:05 AM, Florian Meier <florian.meier at koalo.de>wrote:
> Thank you very much for your review.
>
> On 02/17/2014 02:50 AM, Abhyankar, Shrirang G. wrote:
> > Hi Florian,
> > I was able to run your example with a small change (used PetscMalloc
> > instead of PetscMalloc1).
>
> That is strange. When I use PetscMalloc, it says "Memory corruption!".
Run under valgrind. There is a bug in the code.
> > However, the solution on more than 1 processor
> > is different presumably because of some work "TODO" for the ghost
> vertices.
>
> I do not really understand these ghost vertices. I thought it would be
> ok to just skip those when I iterate over the vertices, but not when I
> iterate over the edges, but apparently that is not the case.
Here is how the parallelism works. We partition the edges into disjoint
sets, one for
each process. Vertices shared by edges on different processes are
duplicated, with
the "ghost" vertex sitting on all processes that do not own the vertex.
> > Answers to some of your questions in the code:
> > I) You can change the non-zero pattern of the Jacobian matrix. But note
> > that if the network topology then I think the circuit would need to be
> > re-built.
>
> I do not know if I really change the non-zero pattern. What I mean is
> the following: Some entries in the Jacobian MIGHT be zero not because of
> the structure of the problem is changing, but just because the partial
> derivative for a special x is (near) zero by chance.
The non-zero pattern of the matrix changes only if your topology changes
(or the
discretization I guess).
> > II) If some of the Jacobian matrix entries are constant then you can use
> > MatStoreValues
> >
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatStoreVal
> > ues.html#MatStoreValues to set the constant terms once before the Newton
> > loop and retrieve these via MatRetrieveValues
> >
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatRetrieve
> > Values.html#MatRetrieveValues during the Jacobian evaluation.
>
> That looks great! Do I have to set the non-constant elements in
> MatStoreValues, too, so that the structure does not change?
I would wait to optimize like this until everything works.
Matt
> > III) The Jacobian matrix for your application has a 2 X 2 structure. So
> > you can set the entires in each block through one MatSetValues call.
>
> It has? I would say that is not the case since all partial derivatives
> might be non-zero when there is a relation (look into the PDF I send
> you). But maybe I just do not know what you mean by 2 X 2 structure or
> what is mean to NOT have a 2 X 2 structure.
>
> > IV) If you want to reuse some term that calculated in the function
> > evaluation (Pfree) during the Jacobian evaluation then you can either use
> > your own user context or save the term in the component array.
>
> That is a good idea. I will try to save it into the component array. Can
> I assume that FormFunction always gets called before FormJacobian with
> the same x?
>
> > Thanks once again for contributing this example. Once you have fixed the
> > ghost vertex part and get the same solution with different number of
> > processors, I'll add it to the PETSc repository. Let us know if you have
> > any more questions.
>
> Thanks again,
> Florian
>
> > On 2/14/14 11:25 AM, "Florian Meier" <florian.meier at koalo.de> wrote:
> >
> >> Ok, that was just a stupid typo.
> >> I have now build the first attempt of an example.
> >> I am sure that can be programmed in a much better way,
> >> so I would be very glad about any comment. What is your
> >> preferred review procedure?
> >>
> >> Greetings,
> >> Florian
> >>
> >> static char help[] = "This example demonstrates the DMCircuit interface
> >> for a radio network flow problem.\n\
> >> The available solver options are in the
> >> alohaoptions file.\n\
> >> The file links.txt contains the network
> >> description.\n\
> >> Run this program: mpiexec -n <n> ./PF\n
> \
> >> mpiexec -n <n> ./PF -pfdata <filename>\n";
> >>
> >> /**
> >> This computes a greatly simplified model of a radio network.
> >> Each node in the network generates packets conforming to a poisson
> >> distribution (given rate R_gen).
> >> These packets are routed along fixed paths to a sink. When two packets
> >> are send at the same time,
> >> they are dropped. No channel sensing or retransmission takes place.
> >>
> >> Each active link (i.e. each radio link that will eventually transmit
> >> data) is modeled as a vertex.
> >> Between those vertices there exist two kinds of relations:
> >> A relation of type REL_INTERFERENCE(affected, source), implies that
> >> whenever the source link tries
> >> to send at the same time than the affected link, the packet will collide
> >> and is not received properly.
> >> A relation of type REL_INFLOW(affected, source) implies that each packet
> >> that is send via the source
> >> link and does not collide, adds up to the amount of packets send via the
> >> affected link.
> >>
> >> The overall rate of accessing a radio link i is
> >> R_access,i = R_gen,i + SUM_OVER_INFLOWS( R_success,j )
> >>
> >> The probability that no other link interfering with i accesses the
> >> channel at a given point in time is
> >> P_free,i = PRODUCT_OVER_INTERFERERS( 1 - R_access,j )
> >>
> >> Finally, the overall rate of successful transmissions over link i is
> >> R_success,i = R_access,i * P_free,i
> >>
> >> The input file is structured as follows:
> >> The first line contains L, the number of links, I, the number of
> >> interferences and F, the number of inflow relations.
> >> Then L lines follow with the packet generation rate for each link.
> >> The subsequent I lines describe the affected and the source link of the
> >> respective interference relation.
> >> The final F lines represent the inflow relations.
> >> */
> >>
> >> #include <petsc.h>
> >> #include <petscdmcircuit.h>
> >>
> >> #include <fstream>
> >> #include <iostream>
> >>
> >> using namespace std;
> >>
> >> enum {
> >> REL_INTERFERENCE,
> >> REL_INFLOW
> >> };
> >>
> >> enum {
> >> VAR_ACCESS,
> >> VAR_SUCCESS,
> >> VAR_NVARS
> >> };
> >>
> >> struct _p_LINKDATA{
> >> PetscScalar packet_generation;
> >> };
> >>
> >> typedef struct _p_LINKDATA *LINKDATA;
> >>
> >> struct _p_RELATIONDATA{
> >> PetscInt source;
> >> PetscInt affected;
> >> PetscInt type;
> >> };
> >>
> >> typedef struct _p_RELATIONDATA *RELATIONDATA;
> >>
> >> typedef struct {
> >> PetscInt nlinks,nrelations; /* # of nodes,relations */
> >> LINKDATA links;
> >> RELATIONDATA relations;
> >> } PFDATA;
> >>
> >>
> >> PetscMPIInt rank;
> >>
> >> #undef __FUNCT__
> >> #define __FUNCT__ "FormFunction"
> >> PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
> >> {
> >> PetscErrorCode ierr;
> >> DM circuitdm;
> >> Vec localX,localF;
> >> PetscInt e;
> >> PetscInt v,vStart,vEnd,vaffected,vsource;
> >> const PetscScalar *xarr;
> >> PetscScalar *farr;
> >> PetscInt offset,offsetsource,offsetlink,offsetrel;
> >> DMCircuitComponentGenericDataType *arr;
> >>
> >> PetscFunctionBegin;
> >> ierr = SNESGetDM(snes,&circuitdm);CHKERRQ(ierr);
> >> ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >> ierr = DMGetLocalVector(circuitdm,&localF);CHKERRQ(ierr);
> >> ierr = VecSet(F,0.0);CHKERRQ(ierr);
> >>
> >> ierr =
> >> DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >> ierr =
> >> DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >>
> >> ierr =
> >> DMGlobalToLocalBegin(circuitdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);
> >> ierr =
> >> DMGlobalToLocalEnd(circuitdm,F,INSERT_VALUES,localF);CHKERRQ(ierr);
> >>
> >> ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
> >> ierr = VecGetArray(localF,&farr);CHKERRQ(ierr);
> >>
> >> ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr);
> >> ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr);
> >>
> >> for (v=vStart; v < vEnd; v++) {
> >> PetscBool ghostvtex;
> >> ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr);
> >> if (ghostvtex) {
> >> continue; // TODO is that ok?
> >> }
> >>
> >> PetscInt keyv;
> >> ierr =
> >>
> DMCircuitGetComponentTypeOffset(circuitdm,v,0,&keyv,&offsetlink);CHKERRQ(i
> >> err);
> >>
> >> LINKDATA link = (LINKDATA)(arr+offsetlink);
> >>
> >> PetscScalar inflow = link->packet_generation;
> >> PetscScalar Pfree = 1;
> >>
> >> PetscInt nconnedges;
> >> const PetscInt *connedges;
> >> ierr =
> >>
> DMCircuitGetSupportingEdges(circuitdm,v,&nconnedges,&connedges);CHKERRQ(ie
> >> rr);
> >> for (PetscInt i = 0; i < nconnedges; i++) {
> >> e = connedges[i];
> >>
> >> const PetscInt *cone;
> >> ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr);
> >> vaffected = cone[0];
> >> vsource = cone[1];
> >>
> >> if (vaffected == v) {
> >> PetscInt keye;
> >> ierr =
> >>
> DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie
> >> rr);
> >> RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel);
> >>
> >> ierr =
> >>
> DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr);
> >>
> >> switch (relation->type) {
> >> case REL_INTERFERENCE:
> >> Pfree *= 1 - xarr[offsetsource+VAR_ACCESS];
> >> break;
> >> case REL_INFLOW:
> >> inflow += xarr[offsetsource+VAR_SUCCESS];
> >> break;
> >> }
> >> }
> >> }
> >>
> >> ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr);
> >> farr[offset+VAR_ACCESS] = inflow - xarr[offset+VAR_ACCESS];
> >> farr[offset+VAR_SUCCESS] = xarr[offset+VAR_ACCESS]*Pfree -
> >> xarr[offset+VAR_SUCCESS];
> >> }
> >>
> >> ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
> >> ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr);
> >> ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >>
> >> ierr =
> >> DMLocalToGlobalBegin(circuitdm,localF,ADD_VALUES,F);CHKERRQ(ierr);
> >> ierr = DMLocalToGlobalEnd(circuitdm,localF,ADD_VALUES,F);CHKERRQ(ierr);
> >> ierr = DMRestoreLocalVector(circuitdm,&localF);CHKERRQ(ierr);
> >>
> >> PetscFunctionReturn(0);
> >> }
> >>
> >>
> >> #undef __FUNCT__
> >> #define __FUNCT__ "FormJacobian"
> >> PetscErrorCode FormJacobian(SNES snes,Vec X, Mat *J,Mat
> >> *Jpre,MatStructure *flg,void *appctx)
> >> {
> >> PetscErrorCode ierr;
> >> DM circuitdm;
> >> Vec localX;
> >> PetscInt e;
> >> PetscInt v,vStart,vEnd,vaffected,vsource;
> >> const PetscScalar *xarr;
> >> PetscInt offsetrel;
> >> DMCircuitComponentGenericDataType *arr;
> >> PetscInt row[1],col[1];
> >> PetscScalar values[1];
> >>
> >> PetscFunctionBegin;
> >> *flg = SAME_NONZERO_PATTERN; // TODO ok for this problem?
> >> ierr = MatZeroEntries(*J);CHKERRQ(ierr);
> >>
> >> ierr = SNESGetDM(snes,&circuitdm);CHKERRQ(ierr);
> >> ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >>
> >> ierr =
> >> DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >> ierr =
> >> DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >>
> >> ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr);
> >>
> >> ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr);
> >> ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr);
> >>
> >> for (v=vStart; v < vEnd; v++) {
> >> PetscInt offset,goffset,offsetsource,goffsetsource;
> >> PetscBool ghostvtex;
> >>
> >> ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr);
> >> if (ghostvtex) {
> >> continue; // TODO is that ok?
> >> }
> >>
> >> ierr =
> >> DMCircuitGetVariableGlobalOffset(circuitdm,v,&goffset);CHKERRQ(ierr);
> >>
> >> // TODO some derivatives are constant, can this be handled in
> >> SetInitialValues?
> >>
> >> // TODO can I combine these two MatSetValues?
> >> row[0] = goffset+VAR_ACCESS;
> >> col[0] = goffset+VAR_ACCESS;
> >> values[0] = -1.0;
> >> ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
> >> row[0] = goffset+VAR_SUCCESS;
> >> col[0] = goffset+VAR_SUCCESS;
> >> values[0] = -1.0;
> >> ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
> >>
> >> // Calculate Pfree
> >> // TODO Pfree was already calculated in FromFunction, can I access
> >> that?
> >> PetscScalar Pfree = 1;
> >>
> >> PetscInt nconnedges;
> >> const PetscInt *connedges;
> >> ierr =
> >>
> DMCircuitGetSupportingEdges(circuitdm,v,&nconnedges,&connedges);CHKERRQ(ie
> >> rr);
> >> for (PetscInt i = 0; i < nconnedges; i++) {
> >> e = connedges[i];
> >>
> >> const PetscInt *cone;
> >> ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr);
> >> vaffected = cone[0];
> >> vsource = cone[1];
> >>
> >> if (vaffected == v) {
> >> PetscInt keye;
> >> ierr =
> >>
> DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie
> >> rr);
> >> RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel);
> >>
> >> ierr =
> >>
> DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr);
> >>
> >> if (relation->type == REL_INTERFERENCE) {
> >> Pfree *= 1 - xarr[offsetsource+VAR_ACCESS];
> >> }
> >> }
> >> }
> >>
> >> row[0] = goffset+VAR_SUCCESS;
> >> col[0] = goffset+VAR_ACCESS;
> >> values[0] = Pfree;
> >> ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
> >>
> >> // Set the derivatives
> >> ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr);
> >> PetscScalar inflow = xarr[offset+VAR_ACCESS];
> >>
> >> for (PetscInt i = 0; i < nconnedges; i++) {
> >> e = connedges[i];
> >>
> >> const PetscInt *cone;
> >> ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr);
> >> vaffected = cone[0];
> >> vsource = cone[1];
> >>
> >> if (vaffected == v) {
> >> PetscInt keye;
> >> ierr =
> >>
> DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie
> >> rr);
> >> RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel);
> >>
> >> ierr =
> >>
> DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr);
> >> ierr =
> >>
> DMCircuitGetVariableGlobalOffset(circuitdm,vsource,&goffsetsource);CHKERRQ
> >> (ierr);
> >>
> >> switch (relation->type) {
> >> case REL_INTERFERENCE:
> >> // TODO Ok to set the entries one by one or is there a better
> >> way?
> >> row[0] = goffset+VAR_SUCCESS;
> >> col[0] = goffsetsource+VAR_ACCESS;
> >> values[0] =
> -inflow*(Pfree/(1-xarr[offsetsource+VAR_ACCESS]));
> >> ierr =
> >> MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
> >> break;
> >> case REL_INFLOW:
> >> row[0] = goffset+VAR_ACCESS;
> >> col[0] = goffsetsource+VAR_SUCCESS;
> >> values[0] = 1.0;
> >> ierr =
> >> MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr);
> >> break;
> >> }
> >> }
> >> }
> >> }
> >>
> >> ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr);
> >> ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >>
> >> ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> >> ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
> >>
> >> PetscFunctionReturn(0);
> >> }
> >>
> >>
> >>
> >> #undef __FUNCT__
> >> #define __FUNCT__ "SetInitialValues"
> >> PetscErrorCode SetInitialValues(DM circuitdm,Vec X,void *appctx)
> >> {
> >> PetscErrorCode ierr;
> >> PetscInt v, vStart, vEnd;
> >> Vec localX;
> >> PetscScalar *xarr;
> >> DMCircuitComponentGenericDataType *arr;
> >>
> >> PetscFunctionBegin;
> >> ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr);
> >>
> >> ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >>
> >> ierr = VecSet(X,0.0);CHKERRQ(ierr);
> >> ierr =
> >> DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >> ierr =
> >> DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
> >>
> >> ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr);
> >> ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr);
> >> for (v = vStart; v < vEnd; v++) {
> >> PetscBool ghostvtex;
> >> ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr);
> >> if (ghostvtex) {
> >> continue; // TODO is that ok?
> >> }
> >>
> >> PetscInt offsetlink, offset;
> >> PetscInt keyv;
> >> ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr);
> >> ierr =
> >>
> DMCircuitGetComponentTypeOffset(circuitdm,v,0,&keyv,&offsetlink);CHKERRQ(i
> >> err);
> >>
> >> LINKDATA link = (LINKDATA)(arr+offsetlink);
> >>
> >> xarr[offset+VAR_ACCESS] = link->packet_generation;
> >> xarr[offset+VAR_SUCCESS] = link->packet_generation;
> >> }
> >>
> >> ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr);
> >> ierr =
> >> DMLocalToGlobalBegin(circuitdm,localX,ADD_VALUES,X);CHKERRQ(ierr);
> >> ierr = DMLocalToGlobalEnd(circuitdm,localX,ADD_VALUES,X);CHKERRQ(ierr);
> >> ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr);
> >> PetscFunctionReturn(0);
> >> }
> >>
> >>
> >> #undef __FUNCT__
> >> #define __FUNCT__ "main"
> >> int main(int argc,char ** argv)
> >> {
> >> string inputFile = "links.txt";
> >> PetscErrorCode ierr;
> >> PFDATA pfdata;
> >> PetscInt numEdges=0,numVertices=0;
> >> int *edges = NULL;
> >> PetscInt i;
> >> DM circuitdm;
> >> PetscInt componentkey[2];
> >> PetscLogStage stage1,stage2;
> >> PetscInt size;
> >>
> >> PetscInitialize(&argc,&argv,"alohaoptions",help);
> >>
> >> ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
> >>
> >> /* Create an empty circuit object */
> >> ierr = DMCircuitCreate(PETSC_COMM_WORLD,&circuitdm);CHKERRQ(ierr);
> >>
> >> /* Register the components in the circuit */
> >> ierr = DMCircuitRegisterComponent(circuitdm,"linkstruct",sizeof(struct
> >> _p_LINKDATA),&componentkey[0]);CHKERRQ(ierr);
> >> ierr =
> >> DMCircuitRegisterComponent(circuitdm,"relationstruct",sizeof(struct
> >> _p_RELATIONDATA),&componentkey[1]);CHKERRQ(ierr);
> >>
> >> ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr);
> >> PetscLogStagePush(stage1);
> >>
> >> /* READ THE DATA */
> >> if (!rank) {
> >> /* Only rank 0 reads the data */
> >> ifstream linksF(inputFile.c_str());
> >>
> >> int ninterferences, ninflows;
> >> linksF >> pfdata.nlinks >> ninterferences >> ninflows;
> >>
> >> numVertices = pfdata.nlinks;
> >> numEdges = pfdata.nrelations = ninterferences + ninflows;
> >>
> >> ierr = PetscMalloc1(pfdata.nlinks,&pfdata.links);CHKERRQ(ierr);
> >> ierr =
> >> PetscMalloc1(pfdata.nrelations,&pfdata.relations);CHKERRQ(ierr);
> >>
> >> for(int i = 0; i < numVertices; i++) {
> >> linksF >> pfdata.links[i].packet_generation;
> >> }
> >>
> >> ierr = PetscMalloc1(2*numEdges,&edges);CHKERRQ(ierr);
> >>
> >> for(int i = 0; i < numEdges; i++) {
> >> linksF >> pfdata.relations[i].affected >>
> >> pfdata.relations[i].source;
> >>
> >> pfdata.relations[i].type = (i < ninterferences) ? REL_INTERFERENCE
> >> : REL_INFLOW;
> >>
> >> edges[2*i] = pfdata.relations[i].affected;
> >> edges[2*i+1] = pfdata.relations[i].source;
> >> }
> >>
> >> linksF.close();
> >> }
> >>
> >> PetscLogStagePop();
> >> MPI_Barrier(PETSC_COMM_WORLD);
> >>
> >> ierr = PetscLogStageRegister("Create circuit",&stage2);CHKERRQ(ierr);
> >> PetscLogStagePush(stage2);
> >> /* Set number of nodes/edges */
> >> ierr =
> >>
> DMCircuitSetSizes(circuitdm,numVertices,numEdges,PETSC_DETERMINE,PETSC_DET
> >> ERMINE);CHKERRQ(ierr);
> >> /* Add edge connectivity */
> >> ierr = DMCircuitSetEdgeList(circuitdm,edges);CHKERRQ(ierr);
> >> /* Set up the circuit layout */
> >> ierr = DMCircuitLayoutSetUp(circuitdm);CHKERRQ(ierr);
> >>
> >> if (!rank) {
> >> ierr = PetscFree(edges);CHKERRQ(ierr);
> >> }
> >>
> >> /* Add circuit components */
> >> PetscInt eStart, eEnd, vStart, vEnd;
> >>
> >> ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr);
> >> for (i = vStart; i < vEnd; i++) {
> >> ierr =
> >>
> DMCircuitAddComponent(circuitdm,i,componentkey[0],&pfdata.links[i-vStart])
> >> ;CHKERRQ(ierr);
> >>
> >> /* Add number of variables */
> >> ierr = DMCircuitAddNumVariables(circuitdm,i,VAR_NVARS);CHKERRQ(ierr);
> >> }
> >>
> >> ierr = DMCircuitGetEdgeRange(circuitdm,&eStart,&eEnd);CHKERRQ(ierr);
> >> for (i = eStart; i < eEnd; i++) {
> >> ierr =
> >>
> DMCircuitAddComponent(circuitdm,i,componentkey[1],&pfdata.relations[i-eSta
> >> rt]);CHKERRQ(ierr);
> >> }
> >>
> >> /* Set up DM for use */
> >> ierr = DMSetUp(circuitdm);CHKERRQ(ierr);
> >>
> >> if (!rank) {
> >> ierr = PetscFree(pfdata.links);CHKERRQ(ierr);
> >> ierr = PetscFree(pfdata.relations);CHKERRQ(ierr);
> >> }
> >>
> >>
> >> ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
> >> if (size > 1) {
> >> DM distcircuitdm;
> >> /* Circuit partitioning and distribution of data */
> >> ierr = DMCircuitDistribute(circuitdm,&distcircuitdm);CHKERRQ(ierr);
> >> ierr = DMDestroy(&circuitdm);CHKERRQ(ierr);
> >> circuitdm = distcircuitdm;
> >> }
> >>
> >> PetscLogStagePop();
> >>
> >> Vec X,F;
> >> ierr = DMCreateGlobalVector(circuitdm,&X);CHKERRQ(ierr);
> >> ierr = VecDuplicate(X,&F);CHKERRQ(ierr);
> >>
> >> Mat J;
> >> ierr = DMCreateMatrix(circuitdm,&J);CHKERRQ(ierr);
> >> ierr =
> >>
> MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
> >>
> >> ierr = SetInitialValues(circuitdm,X,NULL);CHKERRQ(ierr);
> >>
> >> SNES snes;
> >> /* HOOK UP SOLVER */
> >> ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
> >> ierr = SNESSetDM(snes,circuitdm);CHKERRQ(ierr);
> >> ierr = SNESSetFunction(snes,F,FormFunction,NULL);CHKERRQ(ierr);
> >> ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr);
> >> ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
> >>
> >> ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr);
> >> ierr = VecView(X,NULL);CHKERRQ(ierr);
> >>
> >> ierr = VecDestroy(&X);CHKERRQ(ierr);
> >> ierr = VecDestroy(&F);CHKERRQ(ierr);
> >> ierr = MatDestroy(&J);CHKERRQ(ierr);
> >>
> >> ierr = SNESDestroy(&snes);CHKERRQ(ierr);
> >> ierr = DMDestroy(&circuitdm);CHKERRQ(ierr);
> >>
> >> PetscFinalize();
> >> return 0;
> >> }
> >
>
--
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140217/cac69cec/attachment-0001.html>
More information about the petsc-users
mailing list