Hello team, I ran into the address sanitizer error that I hope you could help me with. I don't really know what's wrong with the way the code frees memory. The relevant code file is attached. The line number following domain.cpp specifically referenced to the vector _q, which seems a little odd, since some other vectors are constructed and freed the same way.
==1719==ERROR: AddressSanitizer: attempting free on address which was not malloc()-ed: 0x61f0000076c0 in thread T0 #0 0x7fbf195282ca in __interceptor_free (/usr/lib/x86_64-linux-gnu/libasan.so.2+0x982ca) #1 0x7fbf1706f895 in PetscFreeAlign /home/yyy910805/petsc/src/sys/memory/mal.c:87 #2 0x7fbf1731a898 in VecDestroy_Seq /home/yyy910805/petsc/src/vec/vec/impls/seq/bvec2.c:788 #3 0x7fbf1735f795 in VecDestroy /home/yyy910805/petsc/src/vec/vec/interface/vector.c:408 #4 0x40dd0a in Domain::~Domain() /home/yyy910805/scycle/source/domain.cpp:132 #5 0x40b479 in main /home/yyy910805/scycle/source/main.cpp:242 #6 0x7fbf14d2082f in __libc_start_main (/lib/x86_64-linux-gnu/libc.so.6+0x2082f) #7 0x4075d8 in _start (/home/yyy910805/scycle/source/main+0x4075d8) 0x61f0000076c0 is located 1600 bytes inside of 3220-byte region [0x61f000007080,0x61f000007d14) allocated by thread T0 here: #0 0x7fbf19528b32 in __interceptor_memalign (/usr/lib/x86_64-linux-gnu/libasan.so.2+0x98b32) #1 0x7fbf1706f7e0 in PetscMallocAlign /home/yyy910805/petsc/src/sys/memory/mal.c:41 #2 0x7fbf17073022 in PetscTrMallocDefault /home/yyy910805/petsc/src/sys/memory/mtr.c:183 #3 0x7fbf170710a1 in PetscMallocA /home/yyy910805/petsc/src/sys/memory/mal.c:397 #4 0x7fbf17326fb0 in VecCreate_Seq /home/yyy910805/petsc/src/vec/vec/impls/seq/bvec3.c:35 #5 0x7fbf1736f560 in VecSetType /home/yyy910805/petsc/src/vec/vec/interface/vecreg.c:51 #6 0x7fbf1731afae in VecDuplicate_Seq /home/yyy910805/petsc/src/vec/vec/impls/seq/bvec2.c:807 #7 0x7fbf1735eff7 in VecDuplicate /home/yyy910805/petsc/src/vec/vec/interface/vector.c:379 #8 0x4130de in Domain::setFields() /home/yyy910805/scycle/source/domain.cpp:431 #9 0x40c60a in Domain::Domain(char const*) /home/yyy910805/scycle/source/domain.cpp:57 #10 0x40b433 in main /home/yyy910805/scycle/source/main.cpp:242 #11 0x7fbf14d2082f in __libc_start_main (/lib/x86_64-linux-gnu/libc.so.6+0x2082f) SUMMARY: AddressSanitizer: bad-free ??:0 __interceptor_free ==1719==ABORTING Thanks very much! Yuyun
#include "domain.hpp" #define FILENAME "sbpOps_fc.cpp" using namespace std; // member function definitions including constructor // first type of constructor with 1 parameter Domain::Domain(const char *file) : _file(file),_delim(" = "),_outputDir("data/"), _bulkDeformationType("linearElastic"), _momentumBalanceType("quasidynamic"), _sbpType("mfc_coordTrans"),_operatorType("matrix-based"), _sbpCompatibilityType("fullyCompatible"), _gridSpacingType("variableGridSpacing"), _isMMS(0),_loadICs(0), _inputDir("unspecified_"), _order(4),_Ny(-1),_Nz(-1),_Ly(-1),_Lz(-1), _vL(1e-9), _q(NULL),_r(NULL),_y(NULL),_z(NULL),_y0(NULL),_z0(NULL),_dq(-1),_dr(-1), _bCoordTrans(-1) { #if VERBOSE > 1 std::string funcName = "Domain::Domain(const char *file)"; PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); #endif // load data from file loadData(_file); // check domain size and set grid spacing in y direction if (_Ny > 1) { _dq = 1.0 / (_Ny - 1.0); } else { _dq = 1; } // set grid spacing in z-direction if (_Nz > 1) { _dr = 1.0 / (_Nz - 1.0); } else { _dr = 1; } #if VERBOSE > 2 // each processor prints loaded values to screen PetscMPIInt rank,size; MPI_Comm_size(PETSC_COMM_WORLD,&size); MPI_Comm_rank(PETSC_COMM_WORLD,&rank); for (int Ii = 0; Ii < size; Ii++) { view(Ii); } #endif checkInput(); // perform some basic value checking to prevent NaNs setFields(); setScatters(); #if VERBOSE > 1 PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); #endif } // second type of constructor with 3 parameters Domain::Domain(const char *file,PetscInt Ny, PetscInt Nz) : _file(file),_delim(" = "),_outputDir("data/"), _bulkDeformationType("linearElastic"),_momentumBalanceType("quasidynamic"), _sbpType("mfc_coordTrans"),_operatorType("matrix-based"),_sbpCompatibilityType("fullyCompatible"),_gridSpacingType("variableGridSpacing"), _isMMS(0),_loadICs(0),_inputDir("unspecified_"), _order(4),_Ny(Ny),_Nz(Nz),_Ly(-1),_Lz(-1), _vL(1e-9), _q(NULL),_r(NULL),_y(NULL),_z(NULL),_y0(NULL),_z0(NULL),_dq(-1),_dr(-1), _bCoordTrans(-1) { #if VERBOSE > 1 std::string funcName = "Domain::Domain(const char *file,PetscInt Ny, PetscInt Nz)"; PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); #endif loadData(_file); _Ny = Ny; _Nz = Nz; if (_Ny > 1) { _dq = 1.0/(_Ny-1.0); } else { _dq = 1; } if (_Nz > 1) { _dr = 1.0/(_Nz-1.0); } else { _dr = 1; } #if VERBOSE > 2 // each processor prints loaded values to screen PetscMPIInt rank,size; MPI_Comm_size(PETSC_COMM_WORLD,&size); MPI_Comm_rank(PETSC_COMM_WORLD,&rank); for (int Ii = 0; Ii < size; Ii++) { view(Ii); } #endif checkInput(); // perform some basic value checking to prevent NaNs setFields(); setScatters(); #if VERBOSE > 1 PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); #endif } // destructor Domain::~Domain() { #if VERBOSE > 1 std::string funcName = "Domain::~Domain"; PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); #endif // free memory VecDestroy(&_q); VecDestroy(&_r); VecDestroy(&_y); VecDestroy(&_z); VecDestroy(&_y0); VecDestroy(&_z0); // set map iterator, free memory from VecScatter map<string,VecScatter>::iterator it; for (it = _scatters.begin(); it != _scatters.end(); it++ ) { VecScatterDestroy(&it->second); } #if VERBOSE > 1 PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); #endif } // define loadData function, takes 1 parameter - the filename PetscErrorCode Domain::loadData(const char *file) { PetscErrorCode ierr = 0; PetscMPIInt rank,size; #if VERBOSE > 1 std::string funcName = "Domain::loadData"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif // determines size of the group associated with a communicator MPI_Comm_size(PETSC_COMM_WORLD,&size); // determines rank of the calling processes in the communicator MPI_Comm_rank(PETSC_COMM_WORLD,&rank); // read file inputs ifstream infile(file); string line, var, rhs, rhsFull; size_t pos = 0; while (getline(infile, line)) { istringstream iss(line); pos = line.find(_delim); // find position of the delimiter var = line.substr(0,pos); rhs = ""; if (line.length() > (pos + _delim.length())) { rhs = line.substr(pos+_delim.length(),line.npos); } rhsFull = rhs; // everything after _delim // interpret everything after the appearance of a space on the line as a comment pos = rhs.find(" "); rhs = rhs.substr(0,pos); // set variables, convert string to ints/floats if (var.compare("order") == 0) { _order = atoi(rhs.c_str()); } else if (var.compare("Ny") == 0 && _Ny < 0) { _Ny = atoi(rhs.c_str()); } else if (var.compare("Nz") == 0 && _Nz < 0) { _Nz = atoi(rhs.c_str()); } else if (var.compare("Ly") == 0) { _Ly = atof(rhs.c_str()); } else if (var.compare("Lz") == 0) { _Lz = atof(rhs.c_str()); } else if (var.compare("isMMS") == 0) { _isMMS = 0; std::string temp = rhs; if (temp.compare("yes") == 0 || temp.compare("y") == 0) { _isMMS = 1; } } else if (var.compare("sbpType")==0) { _sbpType = rhs; } else if (var.compare("operatorTYpe")==0) { _operatorType = rhs; } else if (var.compare("sbpCompatibilityType")==0) { _sbpCompatibilityType = rhs; } else if (var.compare("gridSpacingType")==0) { _gridSpacingType = rhs; } else if (var.compare("bulkDeformationType")==0) { _bulkDeformationType = rhs; } else if (var.compare("momentumBalanceType")==0) { _momentumBalanceType = rhs; } else if (var.compare("loadICs")==0) { _loadICs = (int)atof(rhs.c_str()); } else if (var.compare("inputDir")==0) { _inputDir = rhs; } else if (var.compare("bCoordTrans")==0) { _bCoordTrans = atof( rhs.c_str() ); } else if (var.compare("outputDir")==0) { _outputDir = rhs; } else if (var.compare("vL")==0) { _vL = atof( rhs.c_str() ); } } #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif return ierr; } // Specified processor prints scalar/string data members to stdout. PetscErrorCode Domain::view(PetscMPIInt rank) { PetscErrorCode ierr = 0; PetscMPIInt localRank; MPI_Comm_rank(PETSC_COMM_WORLD,&localRank); if (localRank==rank) { #if VERBOSE > 1 std::string funcName = "Domain::view"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif // start printing all the inputs PetscPrintf(PETSC_COMM_SELF,"\n\nrank=%i in Domain::view\n",rank); ierr = PetscPrintf(PETSC_COMM_SELF,"order = %i\n",_order);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"Ny = %i\n",_Ny);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"Nz = %i\n",_Nz);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"Ly = %e\n",_Ly);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"Lz = %e\n",_Lz);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"isMMS = %i\n",_isMMS);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"momBalType = %s\n",_momentumBalanceType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"bulkDeformationType = %s\n",_bulkDeformationType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"sbpType = %s\n",_sbpType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"operatorType = %s\n",_operatorType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"sbpCompatibilityType = %s\n",_sbpCompatibilityType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"gridSpacingType = %s\n",_gridSpacingType.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"outputDir = %s\n",_outputDir.c_str());CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr); #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif } return ierr; } // Check that required fields have been set by the input file PetscErrorCode Domain::checkInput() { PetscErrorCode ierr = 0; #if VERBOSE > 1 std::string funcName = "Domain::checkInput"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif assert(_bulkDeformationType.compare("linearElastic") == 0 || _bulkDeformationType.compare("powerLaw") == 0); assert(_sbpCompatibilityType.compare("fullyCompatible") == 0 || _sbpCompatibilityType.compare("compatible") == 0); if (_bCoordTrans > 0.0) { _gridSpacingType = "variableGridSpacing"; } assert(_gridSpacingType.compare("variableGridSpacing") == 0 || _gridSpacingType.compare("constantGridSpacing") == 0); assert(_momentumBalanceType.compare("quasidynamic") == 0 || _momentumBalanceType.compare("dynamic") == 0 || _momentumBalanceType.compare("quasidynamic_and_dynamic") == 0 || _momentumBalanceType.compare("steadyStateIts") == 0); assert(_order == 2 || _order == 4); assert(_Ly > 0 && _Lz > 0); assert(_dq > 0 && !isnan(_dq)); assert(_dr > 0 && !isnan(_dr)); #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif return ierr; } // Save all scalar fields to text file named domain.txt in output directory. // Also writes out coordinate systems q, r, y, z into respective files in output directory // Note that only the rank 0 processor's values will be saved. PetscErrorCode Domain::write() { PetscErrorCode ierr = 0; #if VERBOSE > 1 std::string funcName = "Domain::write"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif // output scalar fields std::string str = _outputDir + "domain.txt"; PetscViewer viewer; // write into file using PetscViewer PetscViewerCreate(PETSC_COMM_WORLD, &viewer); PetscViewerSetType(viewer, PETSCVIEWERASCII); PetscViewerFileSetMode(viewer, FILE_MODE_WRITE); PetscViewerFileSetName(viewer, str.c_str()); // domain properties ierr = PetscViewerASCIIPrintf(viewer,"order = %i\n",_order);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Ny = %i\n",_Ny);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Nz = %i\n",_Nz);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Ly = %g # (km)\n",_Ly);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"Lz = %g # (km)\n",_Lz);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"isMMS = %i\n",_isMMS);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"momBalType = %s\n",_momentumBalanceType.c_str());CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"bulkDeformationType = %s\n",_bulkDeformationType.c_str());CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"sbpType = %s\n",_sbpType.c_str());CHKERRQ(ierr); // linear solve settings ierr = PetscViewerASCIIPrintf(viewer,"bCoordTrans = %.15e\n",_bCoordTrans);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"outputDir = %s\n",_outputDir.c_str());CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); // get number of processors PetscMPIInt size; MPI_Comm_size(MPI_COMM_WORLD, &size); ierr = PetscViewerASCIIPrintf(viewer,"numProcessors = %i\n",size);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); // output q str = _outputDir + "q"; ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,str.c_str(),FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(_q,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); // output r str = _outputDir + "r"; ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,str.c_str(),FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(_r,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); // output y str = _outputDir + "y"; ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,str.c_str(),FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(_y,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); // output z str = _outputDir + "z"; ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,str.c_str(),FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); ierr = VecView(_z,viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif return ierr; } // construct coordinate transform, setting vectors q, r, y, z PetscErrorCode Domain::setFields() { PetscErrorCode ierr = 0; #if VERBOSE > 1 std::string funcName = "Domain::setFields"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif // generate vector _y with size _Ny*_Nz ierr = VecCreate(PETSC_COMM_WORLD,&_y); CHKERRQ(ierr); ierr = VecSetSizes(_y,PETSC_DECIDE,_Ny*_Nz); CHKERRQ(ierr); ierr = VecSetFromOptions(_y); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) _y, "y"); CHKERRQ(ierr); // duplicate _y into _z, _q, _r VecDuplicate(_y,&_z); PetscObjectSetName((PetscObject) _z, "z"); VecDuplicate(_y,&_q); PetscObjectSetName((PetscObject) _q, "q"); VecDuplicate(_y,&_r); PetscObjectSetName((PetscObject) _r, "r"); // construct coordinate transform PetscInt Ii,Istart,Iend,Jj = 0; PetscScalar *y,*z,*q,*r; ierr = VecGetOwnershipRange(_q,&Istart,&Iend);CHKERRQ(ierr); // return pointers to local data arrays (the processor's portion of vector data) VecGetArray(_y,&y); VecGetArray(_z,&z); VecGetArray(_q,&q); VecGetArray(_r,&r); // set vector entries for q, r (coordinate transform) and y, z (no transform) for (Ii=Istart; Ii<Iend; Ii++) { q[Jj] = _dq*(Ii/_Nz); r[Jj] = _dr*(Ii-_Nz*(Ii/_Nz)); // matrix-based, fully compatible, allows curvilinear coordinate transformation if (_sbpType.compare("mfc_coordTrans") ) { y[Jj] = (_dq*_Ly)*(Ii/_Nz); z[Jj] = (_dr*_Lz)*(Ii-_Nz*(Ii/_Nz)); } else { // hardcoded transformation (not available for z) if (_bCoordTrans > 0) { y[Jj] = _Ly * sinh(_bCoordTrans * q[Jj]) / sinh(_bCoordTrans); } // no transformation y[Jj] = q[Jj]*_Ly; z[Jj] = r[Jj]*_Lz; } Jj++; } // restore arrays VecRestoreArray(_y,&y); VecRestoreArray(_z,&z); VecRestoreArray(_q,&q); VecRestoreArray(_r,&r); // load y and z instead, if provided in input file loadVecFromInputFile(_y,_inputDir,"y"); loadVecFromInputFile(_z,_inputDir,"z"); #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif return ierr; } // scatters values from one vector to another // used to get slip on the fault from the displacement vector, i.e., slip = u(1:Nz); shear stress on the fault from the stress vector sxy; surface displacement; surface heat flux PetscErrorCode Domain::setScatters() { PetscErrorCode ierr = 0; #if VERBOSE > 1 std::string funcName = "Domain::setFields"; ierr = PetscPrintf(PETSC_COMM_WORLD,"Starting %s in %s.\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif // set _y0 to be zero vector with length _Nz VecCreate(PETSC_COMM_WORLD,&_y0); VecSetSizes(_y0,PETSC_DECIDE,_Nz); VecSetFromOptions(_y0); VecSet(_y0,0.0); // set _z0 to be zero vector with length _Ny VecCreate(PETSC_COMM_WORLD,&_z0); VecSetSizes(_z0,PETSC_DECIDE,_Ny); VecSetFromOptions(_z0); VecSet(_z0,0.0); { // set up scatter context to take values for y = 0 from body field and put them on a Vec of size Nz PetscInt *indices; IS is; // index set PetscMalloc1(_Nz,&indices); // we want to scatter from index 0 to _Nz - 1, i.e. take the first _Nz components of the vector to scatter from for (PetscInt Ii = 0; Ii<_Nz; Ii++) { indices[Ii] = Ii; } // creates data structure for an index set containing a list of integers ierr = ISCreateGeneral(PETSC_COMM_WORLD, _Nz, indices, PETSC_COPY_VALUES, &is); // creates vector scatter context, scatters values from _y (at indices is) to _y0 (at indices is) ierr = VecScatterCreate(_y, is, _y0, is, &_scatters["body2L"]); CHKERRQ(ierr); // free memory PetscFree(indices); ISDestroy(&is); } { // set up scatter context to take values for y = Ly from body field and put them on a Vec of size Nz // indices to scatter from PetscInt *fi; IS isf; PetscMalloc1(_Nz,&fi); // we want to scatter from index _Ny*_Nz - _Nz to _Ny*_Nz - 1, i.e. the last _Nz entries of the vector to scatter from for (PetscInt Ii = 0; Ii<_Nz; Ii++) { fi[Ii] = Ii + (_Ny*_Nz-_Nz); } ierr = ISCreateGeneral(PETSC_COMM_WORLD, _Nz, fi, PETSC_COPY_VALUES, &isf); // indices to scatter to PetscInt *ti; IS ist; PetscMalloc1(_Nz,&ti); for (PetscInt Ii = 0; Ii<_Nz; Ii++) { ti[Ii] = Ii; } ierr = ISCreateGeneral(PETSC_COMM_WORLD, _Nz, ti, PETSC_COPY_VALUES, &ist); ierr = VecScatterCreate(_y, isf, _y0, ist, &_scatters["body2R"]); CHKERRQ(ierr); // free memory PetscFree(fi); PetscFree(ti); ISDestroy(&isf); ISDestroy(&ist); } { // set up scatter context to take values for z = 0 from body field and put them on a Vec of size Ny // indices to scatter from IS isf; /* creates a data structure for an index set with a list of evenly spaced integers * locally owned portion of index set has length _Ny * first element of locally owned index set is 0 * change to the next index is _Nz (the stride) * takes indices [0, _Nz, 2*_Nz, ..., (_Ny-1)*_Nz] */ ierr = ISCreateStride(PETSC_COMM_WORLD, _Ny, 0, _Nz, &isf); // indices to scatter to PetscInt *ti; IS ist; PetscMalloc1(_Ny,&ti); // length _Ny for (PetscInt Ii=0; Ii<_Ny; Ii++) { ti[Ii] = Ii; } ierr = ISCreateGeneral(PETSC_COMM_WORLD, _Ny, ti, PETSC_COPY_VALUES, &ist); ierr = VecScatterCreate(_y, isf, _z0, ist, &_scatters["body2T"]); CHKERRQ(ierr); // free memory PetscFree(ti); ISDestroy(&isf); ISDestroy(&ist); } { // set up scatter context to take values for z = Lz from body field and put them on a Vec of size Ny // indices to scatter from IS isf; // takes indices [_Nz - 1, 2*_Nz - 1, ..., _Ny*_Nz - 1] ierr = ISCreateStride(PETSC_COMM_WORLD, _Ny, _Nz - 1, _Nz, &isf); // indices to scatter to PetscInt *ti; IS ist; PetscMalloc1(_Ny,&ti); for (PetscInt Ii = 0; Ii<_Ny; Ii++) { ti[Ii] = Ii; } ierr = ISCreateGeneral(PETSC_COMM_WORLD, _Ny, ti, PETSC_COPY_VALUES, &ist); ierr = VecScatterCreate(_y, isf, _z0, ist, &_scatters["body2B"]); CHKERRQ(ierr); // free memory PetscFree(ti); ISDestroy(&isf); ISDestroy(&ist); } #if VERBOSE > 1 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ending %s in %s\n",funcName.c_str(),FILENAME); CHKERRQ(ierr); #endif return ierr; } // create example vector for testing purposes PetscErrorCode Domain::testScatters() { Vec body; VecDuplicate(_y,&body); PetscInt Istart,Iend,Jj = 0; PetscScalar *bodyA; PetscErrorCode ierr = 0; VecGetOwnershipRange(body,&Istart,&Iend); VecGetArray(body,&bodyA); for (PetscInt Ii = Istart; Ii<Iend; Ii++) { PetscInt Iy = Ii/_Nz; PetscInt Iz = (Ii-_Nz*(Ii/_Nz)); bodyA[Jj] = 10.*Iy + Iz; PetscPrintf(PETSC_COMM_WORLD,"%i %i %g\n",Iy,Iz,bodyA[Jj]); Jj++; } VecRestoreArray(body,&bodyA); // test various mappings // y = 0: mapping to L Vec out; VecDuplicate(_y0,&out); VecScatterBegin(_scatters["body2L"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(_scatters["body2L"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecView(out,PETSC_VIEWER_STDOUT_WORLD); // y = Ly: mapping to R VecDuplicate(_y0,&out); VecSet(out,-1.); VecScatterBegin(_scatters["body2R"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(_scatters["body2R"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecView(out,PETSC_VIEWER_STDOUT_WORLD); // z=0: mapping to T VecDuplicate(_z0,&out); VecSet(out,-1.); VecScatterBegin(_scatters["body2T"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(_scatters["body2T"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecView(out,PETSC_VIEWER_STDOUT_WORLD); // z=Lz: mapping to B VecDuplicate(_z0,&out); VecSet(out,-1.); VecScatterBegin(_scatters["body2B"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecScatterEnd(_scatters["body2B"], body, out, INSERT_VALUES, SCATTER_FORWARD); VecView(out,PETSC_VIEWER_STDOUT_WORLD); // z=Lz: mapping from B to body VecDuplicate(_z0,&out); VecSet(out,-1.); VecScatterBegin(_scatters["body2T"], out, body, INSERT_VALUES, SCATTER_REVERSE); VecScatterEnd(_scatters["body2T"], out, body, INSERT_VALUES, SCATTER_REVERSE); VecView(body,PETSC_VIEWER_STDOUT_WORLD); VecDestroy(&out); VecDestroy(&body); assert(0); return ierr; }