[petsc-users] Modeling a Network with DMPlex?
Abhyankar, Shrirang G.
abhyshr at mcs.anl.gov
Fri Feb 14 13:48:54 CST 2014
Florian,
Thanks for the example! I'll take a look at it and get back to you.
Shri
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;
>}
More information about the petsc-users
mailing list