From 0b62c9de0b4c093b5be1c6bc006205f26ca261e4 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Wed, 16 Aug 2023 10:16:04 +0200 Subject: [PATCH] [EDF27375] : Add InterpKernelDEC.synchronizeWithDefaultValue and retrieveNonFetchedIds methods --- src/ParaMEDMEM/InterpKernelDEC.cxx | 32 ++++++ src/ParaMEDMEM/InterpKernelDEC.hxx | 7 +- src/ParaMEDMEM/InterpolationMatrix.cxx | 76 +++++++++++---- src/ParaMEDMEM/InterpolationMatrix.hxx | 5 + src/ParaMEDMEM/MxN_Mapping.cxx | 27 +++++- src/ParaMEDMEM/MxN_Mapping.hxx | 1 + src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i | 8 ++ src/ParaMEDMEM_Swig/test_InterpKernelDEC.py | 102 ++++++++++++++++++++ 8 files changed, 231 insertions(+), 27 deletions(-) diff --git a/src/ParaMEDMEM/InterpKernelDEC.cxx b/src/ParaMEDMEM/InterpKernelDEC.cxx index ca48b0450..18cf01e7c 100644 --- a/src/ParaMEDMEM/InterpKernelDEC.cxx +++ b/src/ParaMEDMEM/InterpKernelDEC.cxx @@ -165,6 +165,15 @@ namespace MEDCoupling _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. @@ -181,6 +190,29 @@ namespace MEDCoupling } } + MCAuto 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 InterpKernelDEC::retrieveNonFetchedIdsSource() const + { + return _interpolation_matrix->retrieveNonFetchedIdsSource(); + } + + MCAuto 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 diff --git a/src/ParaMEDMEM/InterpKernelDEC.hxx b/src/ParaMEDMEM/InterpKernelDEC.hxx index 653a0e398..eaa87526e 100644 --- a/src/ParaMEDMEM/InterpKernelDEC.hxx +++ b/src/ParaMEDMEM/InterpKernelDEC.hxx @@ -135,13 +135,18 @@ namespace MEDCoupling void release(); void synchronize(); + void synchronizeWithDefaultValue(double val); + MCAuto retrieveNonFetchedIds() const; void recvData(); void recvData(double time); void sendData(); void sendData(double time , double deltatime); void prepareSourceDE() { } void prepareTargetDE() { } - private : + private: + MCAuto retrieveNonFetchedIdsSource() const; + MCAuto retrieveNonFetchedIdsTarget() const; + private: InterpolationMatrix* _interpolation_matrix; }; } diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 91a664add..b419e8d0c 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -849,7 +849,10 @@ namespace MEDCoupling _mapping.prepareSendRecv(); } - + MCAuto InterpolationMatrix::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const + { + return _mapping.retrieveNonFetchedIdsTarget(nbTuples); + } /*! \brief performs t=Ws, where t is the target field, s is the source field @@ -866,12 +869,10 @@ namespace MEDCoupling { mcIdType nbcomp = ToIdType(field.getArray()->getNumberOfComponents()); vector 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 @@ -894,17 +895,42 @@ namespace MEDCoupling if (_target_group.containsMyRank()) { - mcIdType nbelems = field.getArray()->getNumberOfTuples() ; - double* value = const_cast (field.getArray()->getPointer()); - for (mcIdType i=0; ifillWithZero(); } //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 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 InterpolationMatrix::retrieveNonFetchedIdsSource() const + { + MCAuto 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; } @@ -941,17 +967,25 @@ namespace MEDCoupling //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 _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_dft_value); + } } } } diff --git a/src/ParaMEDMEM/InterpolationMatrix.hxx b/src/ParaMEDMEM/InterpolationMatrix.hxx index 5712cccc2..a25cbe7c9 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.hxx +++ b/src/ParaMEDMEM/InterpolationMatrix.hxx @@ -54,11 +54,14 @@ namespace MEDCoupling const mcIdType* distant_elems, const std::string& srcMeth, const std::string& targetMeth); void finishContributionW(ElementLocator& elementLocator); void finishContributionL(ElementLocator& elementLocator); + MCAuto retrieveNonFetchedIdsTarget(mcIdType nbTuples) const; void multiply(MEDCouplingFieldDouble& field) const; + MCAuto 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); @@ -93,6 +96,8 @@ namespace MEDCoupling 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 _row_offsets; std::map, mcIdType > _col_offsets; diff --git a/src/ParaMEDMEM/MxN_Mapping.cxx b/src/ParaMEDMEM/MxN_Mapping.cxx index 5380d74ee..2d2ac386a 100644 --- a/src/ParaMEDMEM/MxN_Mapping.cxx +++ b/src/ParaMEDMEM/MxN_Mapping.cxx @@ -134,6 +134,23 @@ namespace MEDCoupling delete[] recvdispls; } + MCAuto MxN_Mapping::retrieveNonFetchedIdsTarget(mcIdType nbTuples) const + { + MCAuto ret(DataArrayIdType::New()); ret->alloc(0,1); + std::vector 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 @@ -153,7 +170,6 @@ namespace MEDCoupling 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()]; @@ -196,16 +212,17 @@ namespace MEDCoupling //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; icompgetIJ(_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) diff --git a/src/ParaMEDMEM/MxN_Mapping.hxx b/src/ParaMEDMEM/MxN_Mapping.hxx index 3d90c16b9..2f47148a9 100644 --- a/src/ParaMEDMEM/MxN_Mapping.hxx +++ b/src/ParaMEDMEM/MxN_Mapping.hxx @@ -45,6 +45,7 @@ namespace MEDCoupling void addElementFromSource(int distant_proc, mcIdType distant_elem); void prepareSendRecv(); void sendRecv(MEDCouplingFieldDouble& field); + MCAuto retrieveNonFetchedIdsTarget(mcIdType nbTuples) const; void sendRecv(double* sendfield, MEDCouplingFieldDouble& field) const ; void reverseSendRecv(double* recvfield, MEDCouplingFieldDouble& field) const ; diff --git a/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i b/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i index 280b5eed5..624ffadd8 100644 --- a/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i +++ b/src/ParaMEDMEM_Swig/ParaMEDMEMCommon.i @@ -75,6 +75,7 @@ using namespace ICoCo; %newobject MEDCoupling::InterpKernelDEC::_NewWithPG_internal; %newobject MEDCoupling::InterpKernelDEC::_NewWithComm_internal; +%newobject MEDCoupling::InterpKernelDEC::retrieveNonFetchedIds; %newobject MEDCoupling::OverlapDEC::_NewWithComm_internal; %feature("unref") ParaSkyLineArray "$this->decrRef();" @@ -296,6 +297,7 @@ namespace MEDCoupling void release(); void synchronize(); + void synchronizeWithDefaultValue(double val); void recvData(); void recvData(double time); void sendData(); @@ -322,6 +324,12 @@ namespace MEDCoupling { return new InterpKernelDEC(src_ids,trg_ids, *(MPI_Comm*)another_comm); // I know, ugly cast ... } + + DataArrayIdType *retrieveNonFetchedIds() const + { + MCAuto ret = self->retrieveNonFetchedIds(); + return ret.retn(); + } } }; diff --git a/src/ParaMEDMEM_Swig/test_InterpKernelDEC.py b/src/ParaMEDMEM_Swig/test_InterpKernelDEC.py index 3167bf943..4b5b289d7 100755 --- a/src/ParaMEDMEM_Swig/test_InterpKernelDEC.py +++ b/src/ParaMEDMEM_Swig/test_InterpKernelDEC.py @@ -255,6 +255,108 @@ class ParaMEDMEM_IK_DEC_Tests(unittest.TestCase): 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() -- 2.39.2