[petsc-users] Modeling a Network with DMPlex?
Florian Meier
florian.meier at koalo.de
Mon Feb 17 08:05:11 CST 2014
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!".
> 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.
> 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.
> 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?
> 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;
>> }
>
More information about the petsc-users
mailing list