_interpolation_matrix->prepare();
}
+ /*!
+ * Set a default value for non fetched entities
+ */
+ void InterpKernelDEC::synchronizeWithDefaultValue(double val)
+ {
+ this->synchronize();
+ if(_interpolation_matrix )
+ _interpolation_matrix->setDefaultValue(val);
+ }
/*!
Receives the data whether the processor is on the working side or on the lazy side. It must match a \a sendData() call on the other side.
}
}
+ MCAuto<DataArrayIdType> InterpKernelDEC::retrieveNonFetchedIds() const
+ {
+ if( _source_group->containsMyRank() )
+ {
+ return this->retrieveNonFetchedIdsSource();
+ }
+ if( _target_group->containsMyRank() )
+ {
+ return this->retrieveNonFetchedIdsTarget();
+ }
+ THROW_IK_EXCEPTION("Not detected side of rank !");
+ }
+
+ MCAuto<DataArrayIdType> InterpKernelDEC::retrieveNonFetchedIdsSource() const
+ {
+ return _interpolation_matrix->retrieveNonFetchedIdsSource();
+ }
+
+ MCAuto<DataArrayIdType> InterpKernelDEC::retrieveNonFetchedIdsTarget() const
+ {
+ mcIdType nbTuples = _local_field->getField()->getNumberOfTuplesExpected();
+ return _interpolation_matrix->retrieveNonFetchedIdsTarget(nbTuples);
+ }
/*!
Receives the data at time \a time in asynchronous mode. The value of the field
void release();
void synchronize();
+ void synchronizeWithDefaultValue(double val);
+ MCAuto<DataArrayIdType> retrieveNonFetchedIds() const;
void recvData();
void recvData(double time);
void sendData();
void sendData(double time , double deltatime);
void prepareSourceDE() { }
void prepareTargetDE() { }
- private :
+ private:
+ MCAuto<DataArrayIdType> retrieveNonFetchedIdsSource() const;
+ MCAuto<DataArrayIdType> retrieveNonFetchedIdsTarget() const;
+ private:
InterpolationMatrix* _interpolation_matrix;
};
}
_mapping.prepareSendRecv();
}
-
+ MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const
+ {
+ return _mapping.retrieveNonFetchedIdsTarget(nbTuples);
+ }
/*!
\brief performs t=Ws, where t is the target field, s is the source field
{
mcIdType nbcomp = ToIdType(field.getArray()->getNumberOfComponents());
vector<double> target_value(_col_offsets.size()* nbcomp,0.0);
-
//computing the matrix multiply on source side
if (_source_group.containsMyRank())
{
mcIdType nbrows = ToIdType(_coeffs.size());
-
// performing W.S
// W is the intersection matrix
// S is the source vector
if (_target_group.containsMyRank())
{
- mcIdType nbelems = field.getArray()->getNumberOfTuples() ;
- double* value = const_cast<double*> (field.getArray()->getPointer());
- for (mcIdType i=0; i<nbelems*nbcomp; i++)
- {
- value[i]=0.0;
- }
+ field.getArray()->fillWithZero();
}
//on source side : sending T=VT^(-1).(W.S)
//on target side :: receiving T and storing it in field
- _mapping.sendRecv(&target_value[0],field);
+ _mapping.sendRecv(target_value.data(),field);
+
+ if( _target_group.containsMyRank() )
+ {
+ if( this->_presence_dft_value )
+ {
+ const MCAuto<DataArrayIdType> nonFetchedEntities = _mapping.retrieveNonFetchedIdsTarget(field.getArray()->getNumberOfTuples());
+ double *fieldPtr( field.getArray()->getPointerSilent() );
+ for( const mcIdType *eltId = nonFetchedEntities->begin() ; eltId != nonFetchedEntities->end() ; ++eltId)
+ std::fill( fieldPtr + (*eltId)*nbcomp, fieldPtr + ((*eltId)+1)*nbcomp, this->_dft_value );
+ }
+ }
+ }
+
+ MCAuto<DataArrayIdType> InterpolationMatrix::retrieveNonFetchedIdsSource() const
+ {
+ MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+ if (_source_group.containsMyRank())
+ {
+ mcIdType nbrows = ToIdType( _coeffs.size() );
+ for (mcIdType irow = 0; irow < nbrows; irow++)
+ {
+ if( _row_offsets[irow+1] == _row_offsets[irow] )
+ {
+ ret->pushBackSilent( irow );
+ }
+ }
+ }
+ else
+ THROW_IK_EXCEPTION("InterpolationMatrix::retrieveNonFetchedIdsSource : not supposed to be called target side !");
+ return ret;
}
//T is the target vector
for (mcIdType irow = 0; irow < nbrows; irow++)
{
- for (mcIdType icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
- {
- int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
- double value = _coeffs[irow][icol-_row_offsets[irow]].second;
- double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
- for (std::size_t icomp=0; icomp<nbcomp; icomp++)
- {
- double coeff_row = source_value[colid*nbcomp+icomp];
- array[irow*nbcomp+icomp] += value*coeff_row/deno;
- }
- }
+ if( _row_offsets[irow+1] > _row_offsets[irow] )
+ {
+ for (mcIdType icol = _row_offsets[irow]; icol < _row_offsets[irow+1]; icol++)
+ {
+ int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
+ double value = _coeffs[irow][icol-_row_offsets[irow]].second;
+ double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
+ for (std::size_t icomp=0; icomp<nbcomp; icomp++)
+ {
+ double coeff_row = source_value[colid*nbcomp+icomp];
+ array[irow*nbcomp+icomp] += value*coeff_row/deno;
+ }
+ }
+ }
+ else
+ {
+ if( _presence_dft_value )
+ std::fill(array+irow*nbcomp,array+(irow+1)*nbcomp,this->_dft_value);
+ }
}
}
}
const mcIdType* distant_elems, const std::string& srcMeth, const std::string& targetMeth);
void finishContributionW(ElementLocator& elementLocator);
void finishContributionL(ElementLocator& elementLocator);
+ MCAuto<DataArrayIdType> retrieveNonFetchedIdsTarget(mcIdType nbTuples) const;
void multiply(MEDCouplingFieldDouble& field) const;
+ MCAuto<DataArrayIdType> retrieveNonFetchedIdsSource() const;
void transposeMultiply(MEDCouplingFieldDouble& field)const;
void prepare();
mcIdType getNbRows() const { return ToIdType(_row_offsets.size()); }
MPIAccessDEC* getAccessDEC() { return _mapping.getAccessDEC(); }
+ void setDefaultValue(double val) { _presence_dft_value = true; _dft_value = val; }
private:
void computeConservVolDenoW(ElementLocator& elementLocator);
void computeIntegralDenoW(ElementLocator& elementLocator);
private:
bool isSurfaceComputationNeeded(const std::string& method) const;
private:
+ bool _presence_dft_value = false;
+ double _dft_value = 0.0;
const MEDCoupling::ParaFIELD *_source_field;
std::vector<mcIdType> _row_offsets;
std::map<std::pair<int,mcIdType>, mcIdType > _col_offsets;
delete[] recvdispls;
}
+ MCAuto<DataArrayIdType> MxN_Mapping::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const
+ {
+ MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
+ std::vector<bool> hitMachine(nbTuples,false);
+ for (int i=0; i< _recv_proc_offsets[_union_group->size()]; i++)
+ {
+ mcIdType recvId(_recv_ids[i]);
+ hitMachine[recvId] = true;
+ }
+ for( mcIdType ituple = 0 ; ituple < nbTuples ; ++ituple )
+ {
+ if( ! hitMachine[ituple] )
+ ret->pushBackSilent( ituple );
+ }
+ return ret;
+ }
+
/*! Exchanging field data between two groups of processes
*
* \param field MEDCoupling field containing the values to be sent
sendbuf = new double[_sending_ids.size()*nbcomp];
if (_recv_ids.size()>0)
recvbuf = new double[_recv_ids.size()*nbcomp];
-
int* sendcounts = new int[_union_group->size()];
int* senddispls=new int[_union_group->size()];
int* recvcounts=new int[_union_group->size()];
//setting the received values in the field
DataArrayDouble *fieldArr=field.getArray();
- double* recvptr=recvbuf;
+ double* recvptr=recvbuf;
for (int i=0; i< _recv_proc_offsets[_union_group->size()]; i++)
{
+ mcIdType recvId(_recv_ids[i]);
for (int icomp=0; icomp<nbcomp; icomp++)
{
- double temp = fieldArr->getIJ(_recv_ids[i],icomp);
- fieldArr->setIJ(_recv_ids[i],icomp,temp+*recvptr);
+ double temp = fieldArr->getIJ(recvId,icomp);
+ fieldArr->setIJ(recvId,icomp,temp+*recvptr);
recvptr++;
}
- }
+ }
if (sendbuf!=0 && getAllToAllMethod()== Native)
delete[] sendbuf;
if (recvbuf !=0)
void addElementFromSource(int distant_proc, mcIdType distant_elem);
void prepareSendRecv();
void sendRecv(MEDCouplingFieldDouble& field);
+ MCAuto<DataArrayIdType> retrieveNonFetchedIdsTarget(mcIdType nbTuples) const;
void sendRecv(double* sendfield, MEDCouplingFieldDouble& field) const ;
void reverseSendRecv(double* recvfield, MEDCouplingFieldDouble& field) const ;
%newobject MEDCoupling::InterpKernelDEC::_NewWithPG_internal;
%newobject MEDCoupling::InterpKernelDEC::_NewWithComm_internal;
+%newobject MEDCoupling::InterpKernelDEC::retrieveNonFetchedIds;
%newobject MEDCoupling::OverlapDEC::_NewWithComm_internal;
%feature("unref") ParaSkyLineArray "$this->decrRef();"
void release();
void synchronize();
+ void synchronizeWithDefaultValue(double val);
void recvData();
void recvData(double time);
void sendData();
{
return new InterpKernelDEC(src_ids,trg_ids, *(MPI_Comm*)another_comm); // I know, ugly cast ...
}
+
+ DataArrayIdType *retrieveNonFetchedIds() const
+ {
+ MCAuto<DataArrayIdType> ret = self->retrieveNonFetchedIds();
+ return ret.retn();
+ }
}
};
source_group.release()
MPI.COMM_WORLD.Barrier()
+ def test_InterpKernelDEC_default(self):
+ """
+ [EDF27375] : Put a default value when non intersecting case
+ """
+ size = MPI.COMM_WORLD.size
+ rank = MPI.COMM_WORLD.rank
+ if size != 4:
+ print("Should be run on 4 procs!")
+ return
+ nproc_source = 2
+ procs_source = list(range(nproc_source))
+ procs_target = list(range(size - nproc_source, size))
+
+ interface = CommInterface()
+ target_group = MPIProcessorGroup(interface, procs_target)
+ source_group = MPIProcessorGroup(interface, procs_source)
+ dec = InterpKernelDEC(source_group, target_group)
+ #
+ MPI.COMM_WORLD.Barrier()
+ if source_group.containsMyRank():
+ mesh = eval("Source_Proc_{}".format(rank))()
+ nb_local=mesh.getNumberOfCells()
+ field = MEDCouplingFieldDouble(ON_CELLS)
+ field.setNature(IntensiveMaximum)
+ field.setMesh(mesh)
+ arr = DataArrayDouble(nb_local) ; arr.iota() ; arr += rank
+ field.setArray(arr)
+ dec.attachLocalField(field)
+ dec.synchronizeWithDefaultValue(-2000.0)
+ dec.sendData()
+ # target -> source
+ dec.recvData()
+ if rank == 0:
+ self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6,0.6,-2000]),1e-12))
+ self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([2])) )
+ if rank == 1:
+ self.assertTrue(field.getArray().isEqual(DataArrayDouble([1.0,-2000]),1e-12))
+ self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([1])) )
+ else:
+ mesh = eval("Target_Proc_{}".format(rank))()
+ nb_local=mesh.getNumberOfCells()
+ field = MEDCouplingFieldDouble(ON_CELLS)
+ field.setNature(IntensiveMaximum)
+ field.setMesh(mesh)
+ arr = DataArrayDouble(nb_local) ; arr[:] = -20
+ field.setArray(arr)
+ dec.attachLocalField(field)
+ dec.synchronizeWithDefaultValue(-1000.0)
+ dec.recvData()
+ if rank == 2:
+ # matrix S0 / T2 = [[(0,S0,1),(1,S0,1.5)]
+ # IntensiveMaximum => [[(0,S0,1/2.5),(1,S0,1.5/2.5)]
+ #
+ self.assertTrue(field.getArray().isEqual(DataArrayDouble([0.6]),1e-12))
+ self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([])) )
+ if rank == 3:
+ # matrix S1 / T3 = [[],[(0,S1,1.0)],[(0,S1,2.0)],[]]
+ # IntensiveMaximum => [[],[(0,S1,1.0/1.0)],[(0,S1,2.0/2.0)],[]]
+ self.assertTrue(field.getArray().isEqual(DataArrayDouble([-1000.0, 1.0, 1.0, -1000.0]),1e-8))
+ self.assertTrue( dec.retrieveNonFetchedIds().isEqual(DataArrayInt([0,3])) )
+ # target -> source
+ dec.sendData()
+
+ # Some clean up that still needs MPI communication, so to be done before MPI_Finalize()
+ dec.release()
+ target_group.release()
+ source_group.release()
+ MPI.COMM_WORLD.Barrier()
+
+def Source_Proc_0():
+ coo = DataArrayDouble([(0,2),(2,2),(4,2),(0,4),(2,4),(4,4),(0,6),(2,6)])
+ m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_QUAD4,[0,1,4,3])
+ m.insertNextCell(NORM_QUAD4,[1,2,5,4])
+ m.insertNextCell(NORM_QUAD4,[3,4,7,6])
+ return m
+
+def Source_Proc_1():
+ coo = DataArrayDouble([(6,2),(8,2),(10,2),(6,4),(8,4),(10,4)])
+ m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_QUAD4,[0,1,4,3])
+ m.insertNextCell(NORM_QUAD4,[1,2,5,4])
+ return m
+
+def Target_Proc_2():
+ coo = DataArrayDouble([(1,0),(3.5,0),(1,3),(3.5,3)])
+ m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_QUAD4,[0,1,3,2])
+ return m
+
+def Target_Proc_3():
+ coo = DataArrayDouble([(6,0),(7,0),(8,0),(9,0),(10,0),
+ (6,1),(7,1),(9,1),(10,1),
+ (7,3),(8,3),
+ (6,4),(7,4)])
+ m = MEDCouplingUMesh("mesh",2) ; m.setCoords(coo) ; m.allocateCells()
+ m.insertNextCell(NORM_QUAD4,[0,1,6,5])
+ m.insertNextCell(NORM_QUAD4,[1,2,10,9])
+ m.insertNextCell(NORM_QUAD4,[5,6,12,11])
+ m.insertNextCell(NORM_QUAD4,[3,4,8,7])
+ return m
+
if __name__ == "__main__":
unittest.main()
MPI.Finalize()
