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;
}
