diff --git a/applications/solvers/dfSprayFoam/Make/options b/applications/solvers/dfSprayFoam/Make/options
index 0b1293a2..9248840f 100644
--- a/applications/solvers/dfSprayFoam/Make/options
+++ b/applications/solvers/dfSprayFoam/Make/options
@@ -5,6 +5,7 @@ EXE_INC = -std=c++14 \
     -Wno-unused-but-set-variable \
     -Wno-old-style-cast \
     -I. \
+    $(PFLAGS) $(PINC) \
     -I$(FOAM_APP)/solvers/lagrangian/reactingParcelFoam \
     -I$(FOAM_APP)/solvers/compressible/rhoPimpleFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -30,6 +31,7 @@ EXE_INC = -std=c++14 \
     -I$(DF_SRC)/dfCanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
     -I$(DF_SRC)/dfCombustionModels/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
     -I$(CANTERA_ROOT)/include \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \
diff --git a/examples/df0DFoam/zeroD_cubicReactor/H2/pytorchIntegrator/constant/CanteraTorchProperties b/examples/df0DFoam/zeroD_cubicReactor/H2/pytorchIntegrator/constant/CanteraTorchProperties
index ebf63603..dc261039 100644
--- a/examples/df0DFoam/zeroD_cubicReactor/H2/pytorchIntegrator/constant/CanteraTorchProperties
+++ b/examples/df0DFoam/zeroD_cubicReactor/H2/pytorchIntegrator/constant/CanteraTorchProperties
@@ -38,7 +38,7 @@ splittingStrategy    off;
 
 TorchSettings
 {
-    torch            on;
+    torch            off;
     GPU              off;
     log              on;
     torchModel       "HE04_Hydrogen_ESH2_GMS_sub_20221101"; 
diff --git a/src/dfChemistryModel/DNNInferencer/DNNInferencer.cpp b/src/dfChemistryModel/DNNInferencer/DNNInferencer.cpp
index 69aa784e..a008f769 100644
--- a/src/dfChemistryModel/DNNInferencer/DNNInferencer.cpp
+++ b/src/dfChemistryModel/DNNInferencer/DNNInferencer.cpp
@@ -166,7 +166,6 @@ DNNInferencer::DNNInferencer(torch::jit::script::Module torchModel0, torch::jit:
                                313.3717647966624,
                                2.463374792192512e-10},
                               opts);
-    std::cout << "index = " << int(device_.index()) << std::endl;
     std::cout << "load model and parameters successfully" << std::endl;
 }
 
diff --git a/src/dfChemistryModel/Make/options b/src/dfChemistryModel/Make/options
index 80446efa..411866f7 100644
--- a/src/dfChemistryModel/Make/options
+++ b/src/dfChemistryModel/Make/options
@@ -2,6 +2,7 @@
 
 EXE_INC = -std=c++14 \
     -Wno-old-style-cast \
+    $(PFLAGS) $(PINC) \
     $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \
     $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -9,6 +10,7 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
     -I$(DF_SRC)/dfCanteraMixture/lnInclude \
     -I$(CANTERA_ROOT)/include \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
diff --git a/src/dfChemistryModel/dfChemistryModel.C b/src/dfChemistryModel/dfChemistryModel.C
index af5387af..3b086b86 100644
--- a/src/dfChemistryModel/dfChemistryModel.C
+++ b/src/dfChemistryModel/dfChemistryModel.C
@@ -111,6 +111,7 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
         scalar(0.0)
     )
 {
+
 #if defined USE_LIBTORCH || defined USE_PYTORCH
     useDNN = true;
     if (!Qdot_.typeHeaderOk<volScalarField>())
@@ -130,6 +131,7 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     time_getDNNinputs_ = 0;
     time_DNNinference_ = 0;
     time_updateSolutionBuffer_ = 0;
+    time_getProblems_ = 0;
 #endif
 
 #ifdef USE_LIBTORCH
@@ -147,13 +149,9 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     {
         if(!(Pstream::myProcNo() % cores_)) // Now is a master
         {
-
-            Info << "location 0" << endl;
-            Info << "torchModelName1_ = " << torchModelName1_ << endl;
             torch::jit::script::Module torchModel1_ = torch::jit::load(torchModelName1_);
             torch::jit::script::Module torchModel2_ = torch::jit::load(torchModelName2_);
             torch::jit::script::Module torchModel3_ = torch::jit::load(torchModelName3_);
-            Info << "location 1" << endl;
             std::string device_;
             if (gpu_)
             {
@@ -177,6 +175,29 @@ Foam::dfChemistryModel<ThermoType>::dfChemistryModel
     time_python_ = 0;
 #endif
 
+#if defined USE_LIBTORCH || defined USE_PYTORCH
+    // if use torch, create new communicator for solving cvode
+    if (torchSwitch_)
+    {
+        labelList subRank;
+        for (size_t rank = 0; rank < Pstream::nProcs(); rank ++)
+        {
+            if (rank % cores_)
+            {
+                subRank.append(rank);
+            }
+        }
+        cvodeComm = UPstream::allocateCommunicator(UPstream::worldComm, subRank, true);
+        if(Pstream::myProcNo() % cores_)
+        {
+            label sub_rank;
+            MPI_Comm_rank(PstreamGlobals::MPICommunicators_[cvodeComm], &sub_rank);
+            std::cout<<"my ProcessNo in worldComm = " << Pstream::myProcNo() << ' '
+            << "my ProcessNo in cvodeComm = "<<Pstream::myProcNo(cvodeComm)<<std::endl;
+        }
+    }
+#endif
+
     for(const auto& name : CanteraGas_->speciesNames())
     {
         species_.append(name);
@@ -555,6 +576,7 @@ void Foam::dfChemistryModel<ThermoType>::solveSingle
     solution.Qdoti = Qdoti_;
 
     solution.cellid = problem.cellid;
+    solution.local = problem.local;
 }
 
 
@@ -642,7 +664,8 @@ template <class ThermoType>
 Foam::scalar
 Foam::dfChemistryModel<ThermoType>::updateReactionRates
 (
-    const RecvBuffer<ChemistrySolution>& solutions
+    const RecvBuffer<ChemistrySolution>& solutions,
+    DynamicList<ChemistrySolution>& submasterODESolutions
 )
 {
     scalar deltaTMin = great;
@@ -651,14 +674,20 @@ Foam::dfChemistryModel<ThermoType>::updateReactionRates
     {
         for(const auto& solution : array)
         {
+            if (solution.local)
+            {
+                for(label j = 0; j < mixture_.nSpecies(); j++)
+                {
+                    this->RR_[j][solution.cellid] = solution.RRi[j];
+                }
+                this->Qdot_[solution.cellid] = solution.Qdoti;
 
-            for(label j = 0; j < mixture_.nSpecies(); j++)
+                cpuTimes_[solution.cellid] = solution.cpuTime;
+            }
+            else
             {
-                this->RR_[j][solution.cellid] = solution.RRi[j];
+                submasterODESolutions.append(solution);
             }
-            this->Qdot_[solution.cellid] = solution.Qdoti;
-
-            cpuTimes_[solution.cellid] = solution.cpuTime;
         }
     }
 
@@ -758,18 +787,19 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_CVODE
                         << setw(22) << Pstream::myProcNo()
                         << endl;
     }
-
+    DynamicList<ChemistrySolution> List;
     Info<<"=== end solve_CVODE === "<<endl;
-    return updateReactionRates(incomingSolutions);
+    return updateReactionRates(incomingSolutions, List);
 }
 
 #if defined USE_LIBTORCH || defined USE_PYTORCH
 template <class ThermoType>
 template<class DeltaTType>
-Foam::DynamicList<Foam::GpuProblem>
-Foam::dfChemistryModel<ThermoType>::getGPUProblems
+void Foam::dfChemistryModel<ThermoType>::getGPUProblems
 (
-    const DeltaTType& deltaT
+    const DeltaTType &deltaT,
+    Foam::DynamicList<GpuProblem>& GPUproblemList,
+    Foam::DynamicList<ChemistryProblem>& CPUproblemList
 )
 {
     DynamicList<GpuProblem> problemList; //single core TODO:rename it
@@ -795,6 +825,7 @@ Foam::dfChemistryModel<ThermoType>::getGPUProblems
 
         // set problems
         GpuProblem problem(mixture_.nSpecies());
+        ChemistryProblem ode_problem(mixture_.nSpecies());
         problem.cellid = cellI;
         problem.Ti = Ti;
         problem.pi = pi/101325;
@@ -805,31 +836,49 @@ Foam::dfChemistryModel<ThermoType>::getGPUProblems
         problem.rhoi = rhoi;
 
         // choose DNN module
-        if ((Qdot_[cellI] < 3e7) && (T_[cellI] < 2000) && ( T_[cellI] >= 700))//choose1
+        if (((Qdot_[cellI] < 3e7) && (T_[cellI] < 2000) && ( T_[cellI] >= 700)) || (T_[cellI] < 700))//choose1
         {
-            problem.DNNid = 0;
-            problemList.append(problem);
+            // if use CVODE
+            ode_problem.Y = problem.Y;
+            ode_problem.Ti = Ti;
+            ode_problem.pi = pi;
+            ode_problem.rhoi = rhoi;
+            ode_problem.deltaT = deltaT[cellI];
+            ode_problem.cpuTime = cpuTimes_[cellI];
+            ode_problem.cellid = cellI;
+            if (!(Pstream::myProcNo() % cores_)) // submaster
+            {
+                ode_problem.local = false;
+            }
+            CPUproblemList.append(ode_problem);
+            
             selectDNN_[cellI]=0;
             continue;
+
+            // if use DNN
+            // problem.DNNid = 0;
+            // GPUproblemList.append(problem);
+            // continue;
         }
-        if(((Qdot_[cellI] >= 3e7) && (T_[cellI] < 2000)&&(T_[cellI] >= 700))||((Qdot_[cellI] > 7e8) && T_[cellI] > 2000))  //choose2
+        if(((Qdot_[cellI] >= 3e7) && (T_[cellI] < 2000)&&(T_[cellI] >= 700))||((Qdot_[cellI] > 7e8) && T_[cellI] > 2000)) //choose2
         {
             problem.DNNid = 1;
-            problemList.append(problem);
+            GPUproblemList.append(problem);
+            
             selectDNN_[cellI]=1;
             continue;
         }
         if  ((Qdot_[cellI] < 7e8) && (T_[cellI] >= 2000) && (Qdot_[cellI]!=0)) //choose3
         {
             problem.DNNid = 2;
-            problemList.append(problem);
+            GPUproblemList.append(problem);
             selectDNN_[cellI]=2;
             continue;
         }
 
     }
 
-    return problemList;
+    return;
 }
 
 template <class ThermoType>
@@ -1018,7 +1067,16 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
     }
 
     /*=============================gather problems=============================*/
-    DynamicList<GpuProblem> problemList = getGPUProblems(deltaT);
+    std::chrono::steady_clock::time_point start10 = std::chrono::steady_clock::now();
+    DynamicList<GpuProblem> GPUproblemList; //single core TODO:rename it
+    DynamicList<ChemistryProblem> CPUproblemList;
+    getGPUProblems(deltaT, GPUproblemList, CPUproblemList);
+    label flag_mpi_init;
+    MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+    std::chrono::steady_clock::time_point stop10 = std::chrono::steady_clock::now();
+    std::chrono::duration<double> processingTime10 = std::chrono::duration_cast<std::chrono::duration<double>>(stop10 - start10);
+    time_getProblems_ += processingTime10.count();
 
     /*==============================send problems==============================*/
     std::chrono::steady_clock::time_point start2 = std::chrono::steady_clock::now();
@@ -1027,10 +1085,28 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
     if (Pstream::myProcNo() % cores_) //for slave
     {
         UOPstream send((Pstream::myProcNo()/cores_)*cores_, pBufs);// sending problem to master
-        send << problemList;
+        send << GPUproblemList;
     }
     pBufs.finishedSends();
 
+    /*==============================send CVODE problems from submaster to neighbour==============================*/
+    PstreamBuffers pBufs1(Pstream::commsTypes::nonBlocking);
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        UOPstream send((Pstream::myProcNo() + 1), pBufs1);// sending CPUproblems to neighbour
+        send << CPUproblemList;
+    }
+    pBufs1.finishedSends();
+    if ((Pstream::myProcNo() % cores_) == 1) // neighbour of submaster
+    {
+        DynamicList<ChemistryProblem> CPUproblemList_submaster;
+        UIPstream recv((Pstream::myProcNo() - 1), pBufs1);
+        recv >> CPUproblemList_submaster;
+        CPUproblemList.append(CPUproblemList_submaster);
+    }
+
+    /*========================================================================================================*/
+
     DynamicBuffer<GpuSolution> solutionBuffer;
 
     std::chrono::steady_clock::time_point stop2 = std::chrono::steady_clock::now();
@@ -1048,7 +1124,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         DynamicBuffer<GpuProblem> problemBuffer(cores_);//each submaster init a local problemBuffer TODO:rename it
 
         /*==============================gather problems==============================*/
-        problemBuffer[0] = problemList; //problemList of submaster get index 0
+        problemBuffer[0] = GPUproblemList; //problemList of submaster get index 0
         problemSize += problemBuffer[0].size();
 
         for (label i = 1; i < cores_; i++)
@@ -1077,7 +1153,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         getDNNinputs(problemBuffer, outputLength, DNNinputs, cellIDBuffer, problemCounter);
         std::chrono::steady_clock::time_point stop5 = std::chrono::steady_clock::now();
         std::chrono::duration<double> processingTime5 = std::chrono::duration_cast<std::chrono::duration<double>>(stop5 - start5);
-        // std::cout << "getDNNinputsTime = " << processingTime5.count() << std::endl;
+        std::cout << "getDNNinputsTime = " << processingTime5.count() << std::endl;
         time_getDNNinputs_ += processingTime5.count();
 
         /*=============================inference via DNNInferencer=============================*/
@@ -1087,7 +1163,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
 
         std::chrono::steady_clock::time_point stop7 = std::chrono::steady_clock::now();
         std::chrono::duration<double> processingTime7 = std::chrono::duration_cast<std::chrono::duration<double>>(stop7 - start7);
-        // std::cout << "DNNinferenceTime = " << processingTime7.count() << std::endl;
+        std::cout << "DNNinferenceTime = " << processingTime7.count() << std::endl;
         time_DNNinference_ += processingTime7.count();
 
         /*=============================construct solutions=============================*/
@@ -1097,7 +1173,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
 
         std::chrono::steady_clock::time_point stop6 = std::chrono::steady_clock::now();
         std::chrono::duration<double> processingTime6 = std::chrono::duration_cast<std::chrono::duration<double>>(stop6 - start6);
-        // std::cout << "updateSolutionBufferTime = " << processingTime6.count() << std::endl;
+        std::cout << "updateSolutionBufferTime = " << processingTime6.count() << std::endl;
         time_updateSolutionBuffer_ += processingTime6.count();
 
         std::chrono::steady_clock::time_point stop1 = std::chrono::steady_clock::now();
@@ -1106,12 +1182,46 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         time_submaster_ += processingTime1.count();
     }
 
-    /*=============================send and recv solutions=============================*/
+    /*=============================calculates RR with CVODE use DLB=============================*/
+    DynamicList<ChemistrySolution> CPUSolutionList;
+    if (Pstream::myProcNo() % cores_) //for slave
+    {
+        std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
+        DynamicBuffer<ChemistrySolution> incomingSolutions;
+        balancer_.updateState(CPUproblemList, cvodeComm);
+        auto guestProblems = balancer_.balance(CPUproblemList, cvodeComm);
+        auto ownProblems = balancer_.getRemaining(CPUproblemList, cvodeComm);
+        auto ownSolutions = solveList(ownProblems);
+        auto guestSolutions = solveBuffer(guestProblems);
+        incomingSolutions = balancer_.unbalance(guestSolutions, cvodeComm);
+        incomingSolutions.append(ownSolutions);
+        updateReactionRates(incomingSolutions, CPUSolutionList);
+        std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
+        std::chrono::duration<double> processingTime = std::chrono::duration_cast<std::chrono::duration<double>>(stop - start);
+        std::cout << "slaveTime = " << processingTime.count() << std::endl;
+    }
+
+    /*=============================send CPUSolutionList back to submaster=============================*/
+    PstreamBuffers pBufs3(Pstream::commsTypes::nonBlocking);
+
+    if ((Pstream::myProcNo() % cores_) == 1) // neighbour of submaster
+    {
+        UOPstream send((Pstream::myProcNo() - 1), pBufs3);
+        send << CPUSolutionList;
+    }
+    pBufs3.finishedSends();
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        UIPstream recv((Pstream::myProcNo() + 1), pBufs3);// resv CPUproblems from neighbour
+        recv >> CPUSolutionList;
+    }
+
+    /*=============================send and recv GPUSolutions=============================*/
     std::chrono::steady_clock::time_point start4 = std::chrono::steady_clock::now();
 
     DynamicList<GpuSolution> finalList;
     PstreamBuffers pBufs2(Pstream::commsTypes::nonBlocking);
-    if (!(Pstream::myProcNo() % cores_))
+    if (!(Pstream::myProcNo() % cores_)) // submaster
     {
         finalList = solutionBuffer[0];
         for (label i = 1; i < cores_; i++)
@@ -1123,7 +1233,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
     }
     pBufs2.finishedSends();
 
-    if (Pstream::myProcNo() % cores_)
+    if (Pstream::myProcNo() % cores_) // slavers
     {
         UIPstream recv((Pstream::myProcNo()/cores_)*cores_, pBufs2);
         recv >> finalList;
@@ -1145,6 +1255,19 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         }
     }
 
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        for (int cellI = 0; cellI < CPUSolutionList.size(); cellI++)
+        {
+            for (int speciID = 0; speciID < mixture_.nSpecies(); speciID++)
+            {
+                RR_[speciID][CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].RRi[speciID];
+            }
+            Qdot_[CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].Qdoti;
+            cpuTimes_[CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].cpuTime;
+        }
+    }
+
     Info << "=== end solve_DNN === " << endl;
 
     std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
@@ -1326,7 +1449,16 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
 
 
     /*=============================gather problems=============================*/
-    DynamicList<GpuProblem> problemList = getGPUProblems(deltaT);
+    std::chrono::steady_clock::time_point start10 = std::chrono::steady_clock::now();
+    DynamicList<GpuProblem> GPUproblemList; //single core TODO:rename it
+    DynamicList<ChemistryProblem> CPUproblemList;
+    getGPUProblems(deltaT, GPUproblemList, CPUproblemList);
+    label flag_mpi_init;
+    MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+    std::chrono::steady_clock::time_point stop10 = std::chrono::steady_clock::now();
+    std::chrono::duration<double> processingTime10 = std::chrono::duration_cast<std::chrono::duration<double>>(stop10 - start10);
+    time_getProblems_ += processingTime10.count();
 
     /*==============================send problems==============================*/
     std::chrono::steady_clock::time_point start2 = std::chrono::steady_clock::now();
@@ -1335,10 +1467,28 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
     if (Pstream::myProcNo() % cores_) //for slave
     {
         UOPstream send((Pstream::myProcNo()/cores_)*cores_, pBufs);// sending problem to master
-        send << problemList;
+        send << GPUproblemList;
     }
     pBufs.finishedSends();
 
+    /*==============================send CVODE problems from submaster to neighbour==============================*/
+    PstreamBuffers pBufs1(Pstream::commsTypes::nonBlocking);
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        UOPstream send((Pstream::myProcNo() + 1), pBufs1);// sending CPUproblems to neighbour
+        send << CPUproblemList;
+    }
+    pBufs1.finishedSends();
+    if ((Pstream::myProcNo() % cores_) == 1) // neighbour of submaster
+    {
+        DynamicList<ChemistryProblem> CPUproblemList_submaster;
+        UIPstream recv((Pstream::myProcNo() - 1), pBufs1);
+        recv >> CPUproblemList_submaster;
+        CPUproblemList.append(CPUproblemList_submaster);
+    }
+
+    /*========================================================================================================*/
+
     DynamicBuffer<GpuSolution> solutionBuffer;
 
     std::chrono::steady_clock::time_point stop2 = std::chrono::steady_clock::now();
@@ -1356,7 +1506,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         DynamicBuffer<GpuProblem> problemBuffer(cores_);//each submaster init a local problemBuffer TODO:rename it
 
         /*==============================gather problems==============================*/
-        problemBuffer[0] = problemList; //problemList of submaster get index 0
+        problemBuffer[0] = GPUproblemList; //problemList of submaster get index 0
         problemSize += problemBuffer[0].size();
 
         for (label i = 1; i < cores_; i++)
@@ -1436,12 +1586,46 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         time_submaster_ += processingTime1.count();
     }
 
-    /*=============================send and recv solutions=============================*/
+    /*=============================calculates RR with CVODE use DLB=============================*/
+    DynamicList<ChemistrySolution> CPUSolutionList;
+    if (Pstream::myProcNo() % cores_) //for slave
+    {
+        std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
+        DynamicBuffer<ChemistrySolution> incomingSolutions;
+        balancer_.updateState(CPUproblemList, cvodeComm);
+        auto guestProblems = balancer_.balance(CPUproblemList, cvodeComm);
+        auto ownProblems = balancer_.getRemaining(CPUproblemList, cvodeComm);
+        auto ownSolutions = solveList(ownProblems);
+        auto guestSolutions = solveBuffer(guestProblems);
+        incomingSolutions = balancer_.unbalance(guestSolutions, cvodeComm);
+        incomingSolutions.append(ownSolutions);
+        updateReactionRates(incomingSolutions, CPUSolutionList);
+        std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
+        std::chrono::duration<double> processingTime = std::chrono::duration_cast<std::chrono::duration<double>>(stop - start);
+        std::cout << "slaveTime = " << processingTime.count() << std::endl;
+    }
+
+    /*=============================send CPUSolutionList back to submaster=============================*/
+    PstreamBuffers pBufs3(Pstream::commsTypes::nonBlocking);
+
+    if ((Pstream::myProcNo() % cores_) == 1) // neighbour of submaster
+    {
+        UOPstream send((Pstream::myProcNo() - 1), pBufs3);
+        send << CPUSolutionList;
+    }
+    pBufs3.finishedSends();
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        UIPstream recv((Pstream::myProcNo() + 1), pBufs3);// resv CPUproblems from neighbour
+        recv >> CPUSolutionList;
+    }
+
+    /*=============================send and recv GPUSolutions=============================*/
     std::chrono::steady_clock::time_point start4 = std::chrono::steady_clock::now();
 
     DynamicList<GpuSolution> finalList;
     PstreamBuffers pBufs2(Pstream::commsTypes::nonBlocking);
-    if (!(Pstream::myProcNo() % cores_))
+    if (!(Pstream::myProcNo() % cores_)) // submaster
     {
         finalList = solutionBuffer[0];
         for (label i = 1; i < cores_; i++)
@@ -1451,7 +1635,7 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         }
     }
     pBufs2.finishedSends();
-    if (Pstream::myProcNo() % cores_)
+    if (Pstream::myProcNo() % cores_) // slavers
     {
         UIPstream recv((Pstream::myProcNo()/cores_)*cores_, pBufs2);
         recv >> finalList;
@@ -1473,6 +1657,19 @@ Foam::scalar Foam::dfChemistryModel<ThermoType>::solve_DNN(
         }
     }
 
+    if (!(Pstream::myProcNo() % cores_)) // submaster
+    {
+        for (int cellI = 0; cellI < CPUSolutionList.size(); cellI++)
+        {
+            for (int speciID = 0; speciID < mixture_.nSpecies(); speciID++)
+            {
+                RR_[speciID][CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].RRi[speciID];
+            }
+            Qdot_[CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].Qdoti;
+            cpuTimes_[CPUSolutionList[cellI].cellid] = CPUSolutionList[cellI].cpuTime;
+        }
+    }
+
     Info << "=== end solve_DNN === " << endl;
 
     std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
diff --git a/src/dfChemistryModel/dfChemistryModel.H b/src/dfChemistryModel/dfChemistryModel.H
index 618c6fea..ad2a4e55 100644
--- a/src/dfChemistryModel/dfChemistryModel.H
+++ b/src/dfChemistryModel/dfChemistryModel.H
@@ -71,6 +71,7 @@ SourceFiles
 #include "LoadBalancer.H"
 #include "OFstream.H"
 #include "IOmanip.H"
+#include "PstreamGlobals.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -140,6 +141,7 @@ public IOdictionary
     scalar Qdotact2_;
     scalar Tact3_;
     scalar Qdotact3_;
+    label cvodeComm;
 
     // profiling 
     double time_allsolve_;
@@ -150,6 +152,7 @@ public IOdictionary
     double time_getDNNinputs_;
     double time_DNNinference_;
     double time_updateSolutionBuffer_;
+    double time_getProblems_;
 #endif
 
 #ifdef USE_LIBTORCH
@@ -218,7 +221,7 @@ public IOdictionary
         solveBuffer(RecvBuffer<ChemistryProblem>& problems);
 
         //- Update the reaction rates and heat release rates from a list of solutions
-        scalar updateReactionRates(const RecvBuffer<ChemistrySolution>& solutions);
+        scalar updateReactionRates(const RecvBuffer<ChemistrySolution>& solutions, DynamicList<ChemistrySolution>& submasterODESolutions);
 
         //- Create a load balancer object
         LoadBalancer createBalancer();
@@ -239,7 +242,8 @@ public IOdictionary
         /*============= private methods in solve_DNN =============*/
         //- get a list of GPU problems to be solved
         template<class DeltaTType>
-        DynamicList<GpuProblem> getGPUProblems(const DeltaTType& deltaT);
+        void getGPUProblems(const DeltaTType& deltaT, Foam::DynamicList<GpuProblem>& GPUproblemList,
+        Foam::DynamicList<ChemistryProblem>& CPUproblemList);
 
         //- get the input for DNN inference
         void getDNNinputs(const DynamicBuffer<GpuProblem>& problemBuffer, std::vector<label>& outputlength,
@@ -322,6 +326,7 @@ public:
     double time_allsolve() {return time_allsolve_;}
     double time_submaster() {return time_submaster_;}
     double time_sendProblem() {return time_sendProblem_;}
+    double time_getProblem() {return time_getProblems_;}
     double time_RecvProblem() {return time_RecvProblem_;}
     double time_sendRecvSolution() {return time_sendRecvSolution_;}
     double time_getDNNinputs() {return time_getDNNinputs_;}
diff --git a/src/dfChemistryModel/loadBalancing/ChemistryProblem.H b/src/dfChemistryModel/loadBalancing/ChemistryProblem.H
index 123a71ef..4ba67857 100644
--- a/src/dfChemistryModel/loadBalancing/ChemistryProblem.H
+++ b/src/dfChemistryModel/loadBalancing/ChemistryProblem.H
@@ -47,7 +47,7 @@ struct ChemistryProblem
 
     ChemistryProblem() = default;
     ChemistryProblem(label nSpecie)
-        : Y(nSpecie), Ti(0), pi(0), rhoi(0), deltaT(0), cellid(0)
+        : Y(nSpecie), Ti(0), pi(0), rhoi(0), deltaT(0), cellid(0), local(true)
     {
     }
 
@@ -58,6 +58,7 @@ struct ChemistryProblem
     scalar deltaT;
     scalar cpuTime;
     label cellid;
+    bool local;
 
     // TODO: implement!
     bool operator==(const ChemistryProblem& rhs) const
@@ -82,6 +83,7 @@ static inline Ostream& operator<<(Ostream& os, const ChemistryProblem& p)
     os << p.deltaT;
     os << p.cpuTime;
     os << p.cellid;
+    os << p.local;
 
     return os;
 }
@@ -97,6 +99,7 @@ static inline Istream& operator>>(Istream& is, ChemistryProblem& p)
     is >> p.deltaT;
     is >> p.cpuTime;
     is >> p.cellid;
+    is >> p.local;
 
     return is;
 }
diff --git a/src/dfChemistryModel/loadBalancing/ChemistrySolution.H b/src/dfChemistryModel/loadBalancing/ChemistrySolution.H
index 84fc89d7..0527d804 100644
--- a/src/dfChemistryModel/loadBalancing/ChemistrySolution.H
+++ b/src/dfChemistryModel/loadBalancing/ChemistrySolution.H
@@ -48,7 +48,7 @@ struct ChemistrySolution
     ChemistrySolution() = default;
 
     ChemistrySolution(label nspecie)
-        : RRi(nspecie, 0.0), cpuTime(0.0), cellid(0), Qdoti(0.0)
+        : RRi(nspecie, 0.0), cpuTime(0.0), cellid(0), Qdoti(0.0), local(true)
     {
     }
 
@@ -66,6 +66,7 @@ struct ChemistrySolution
     scalar cpuTime;
     label cellid;
     scalar Qdoti;
+    bool local;
 };
 
 //- Serialization for send
@@ -75,6 +76,7 @@ static inline Ostream& operator<<(Ostream& os, const ChemistrySolution& s)
     os << s.cpuTime;
     os << s.cellid;
     os << s.Qdoti;
+    os << s.local;
     return os;
 }
 
@@ -85,6 +87,7 @@ static inline Istream& operator>>(Istream& is, ChemistrySolution& s)
     is >> s.cpuTime;
     is >> s.cellid;
     is >> s.Qdoti;
+    is >> s.local;
     return is;
 }
 
diff --git a/src/dfChemistryModel/loadBalancing/LoadBalancer.C b/src/dfChemistryModel/loadBalancing/LoadBalancer.C
index a4065534..d7a5d094 100644
--- a/src/dfChemistryModel/loadBalancing/LoadBalancer.C
+++ b/src/dfChemistryModel/loadBalancing/LoadBalancer.C
@@ -28,10 +28,11 @@ License
 
 void
 Foam::LoadBalancer::updateState(
-    const DynamicList<ChemistryProblem>& problems)
+    const DynamicList<ChemistryProblem>& problems, const label comm)
 {
-    auto myLoad = computeLoad(problems);
-    auto allLoads = allGather(myLoad);
+    if (Pstream::myProcNo(comm) == -1) return;
+    auto myLoad = computeLoad(problems, comm);
+    auto allLoads = allGather(myLoad, comm);
     auto operations = getOperations(allLoads, myLoad);
     auto info = operationsToInfo(operations, problems, myLoad);
 
diff --git a/src/dfChemistryModel/loadBalancing/LoadBalancer.H b/src/dfChemistryModel/loadBalancing/LoadBalancer.H
index 3bee83b9..90fceb76 100644
--- a/src/dfChemistryModel/loadBalancing/LoadBalancer.H
+++ b/src/dfChemistryModel/loadBalancing/LoadBalancer.H
@@ -79,7 +79,8 @@ public:
     virtual ~LoadBalancer() = default;
 
     //- Given a list of problems, update the balancer state member
-    virtual void updateState(const DynamicList<ChemistryProblem>& problems);
+    virtual void updateState(const DynamicList<ChemistryProblem>& problems,
+    const label comm = UPstream::worldComm);
 
     //- Is load balancing active?
     bool active() const
diff --git a/src/dfChemistryModel/loadBalancing/LoadBalancerBase.C b/src/dfChemistryModel/loadBalancing/LoadBalancerBase.C
index d5ba22fd..0d7413f3 100644
--- a/src/dfChemistryModel/loadBalancing/LoadBalancerBase.C
+++ b/src/dfChemistryModel/loadBalancing/LoadBalancerBase.C
@@ -40,14 +40,16 @@ void Foam::LoadBalancerBase::setState(const BalancerState& state)
 }
 
 Foam::ChemistryLoad
-Foam::LoadBalancerBase::computeLoad(const DynamicList<ChemistryProblem>& problems)
+Foam::LoadBalancerBase::computeLoad(const DynamicList<ChemistryProblem>& problems, 
+const label comm)
 {
     auto lambda = [](scalar sum, const ChemistryProblem& rhs) {
         return sum + rhs.cpuTime;
     };
     scalar sum =
         std::accumulate(problems.begin(), problems.end(), scalar(0), lambda);
-    return ChemistryLoad(Pstream::myProcNo(), sum);
+    
+    return ChemistryLoad(Pstream::myProcNo(comm), sum);
 }
 
 Foam::scalar Foam::LoadBalancerBase::getMean(const DynamicList<ChemistryLoad>& loads)
diff --git a/src/dfChemistryModel/loadBalancing/LoadBalancerBase.H b/src/dfChemistryModel/loadBalancing/LoadBalancerBase.H
index 86af58d5..f857bcf3 100644
--- a/src/dfChemistryModel/loadBalancing/LoadBalancerBase.H
+++ b/src/dfChemistryModel/loadBalancing/LoadBalancerBase.H
@@ -76,7 +76,8 @@ public:
     virtual bool active() const = 0;
 
     //- The load balancing algorithm which each derived class must define
-    virtual void updateState(const DynamicList<ChemistryProblem>& problems) = 0;
+    virtual void updateState(const DynamicList<ChemistryProblem>& problems, 
+    const label comm = UPstream::worldComm) = 0;
 
     //- Calculate the mean of the given list of loads
     static double getMean(const DynamicList<ChemistryLoad>& loads);
@@ -89,11 +90,13 @@ public:
 
     //- Compute the load based on the given list of chemistry problems
     static ChemistryLoad
-    computeLoad(const DynamicList<ChemistryProblem>& problems);
+    computeLoad(const DynamicList<ChemistryProblem>& problems,
+    const label comm = UPstream::worldComm);
 
     //- Gather the data from all ranks
     template <class T>
-    static DynamicList<T> allGather(const T& myData);
+    static DynamicList<T> allGather(const T& myData, 
+    const label comm = UPstream::worldComm);
 
     //- Set the current state of the rank
     void setState(const BalancerState& state);
@@ -106,11 +109,13 @@ public:
 
     //- Given a list of values, split them evenly between MPI processes
     template <class T>
-    RecvBuffer<T> balance(const DynamicList<T>& values) const;
+    RecvBuffer<T> balance(const DynamicList<T>& values, 
+    const label comm = UPstream::worldComm) const;
 
     //- Given a buffer of values, send the values back to their owner ranks
     template <class T>
-    RecvBuffer<T> unbalance(const RecvBuffer<T>& values) const;
+    RecvBuffer<T> unbalance(const RecvBuffer<T>& values,
+    const label comm = UPstream::worldComm) const;
 
     //- Print the current state information
     void printState() const;
@@ -132,48 +137,55 @@ public:
     static RecvBuffer<ET> sendRecv(
         const Indexable&          send_buffer,
         const std::vector<label>& sources,
-        const std::vector<label>& destinations);
+        const std::vector<label>& destinations,
+        const label comm = UPstream::worldComm);
 
     //- Slice an nRemaining size portion from the _end_ of the values
     template <class T>
-    SubList<T> getRemaining(const DynamicList<T>& values)
+    SubList<T> getRemaining(const DynamicList<T>& values, const label comm = UPstream::worldComm)
     {
-
+        if (Pstream::myProcNo(comm) == -1) 
+        {
+            SubList<T> res(values, state_.nRemaining, values.size() - state_.nRemaining);
+            return res.null();
+        }
         return SubList<T>(
             values, state_.nRemaining, values.size() - state_.nRemaining);
     }
 };
 
 template <class T>
-DynamicList<T> LoadBalancerBase::allGather(const T& myData)
+DynamicList<T> LoadBalancerBase::allGather(const T& myData, 
+const label comm)
 {
 
-    DynamicList<T> ret(Pstream::nProcs(), myData);
-    ret[Pstream::myProcNo()] = myData;
+    DynamicList<T> ret(Pstream::nProcs(comm), myData);
+    ret[Pstream::myProcNo(comm)] = myData;
 
     label tag = 1;
-    Pstream::gatherList(ret, tag);
-    Pstream::scatterList(ret, tag);
+    Pstream::gatherList(ret, tag, comm);
+    Pstream::scatterList(ret, tag, comm);
 
     return ret;
 }
 
 template <class T>
-RecvBuffer<T> LoadBalancerBase::balance(const DynamicList<T>& values) const
+RecvBuffer<T> LoadBalancerBase::balance(const DynamicList<T>& values, const label comm) const
 {
 
     return sendRecv<T, SendBuffer<T>>(
         SendBuffer<T>(values, state_.nProblems),
         state_.sources,
-        state_.destinations);
+        state_.destinations,
+        comm);
 }
 
 template <class T>
-RecvBuffer<T> LoadBalancerBase::unbalance(const RecvBuffer<T>& values) const
+RecvBuffer<T> LoadBalancerBase::unbalance(const RecvBuffer<T>& values, const label comm) const
 {
 
     return sendRecv<T, RecvBuffer<T>>(
-        values, state_.destinations, state_.sources);
+        values, state_.destinations, state_.sources, comm);
 }
 
 template <class T>
@@ -194,7 +206,8 @@ template <class ET, class Indexable>
 RecvBuffer<ET> LoadBalancerBase::sendRecv(
     const Indexable&          send_buffer,
     const std::vector<label>& sources,
-    const std::vector<label>& destinations)
+    const std::vector<label>& destinations,
+    const label comm)
 {
 
     RecvBuffer<ET> ret;
@@ -202,21 +215,39 @@ RecvBuffer<ET> LoadBalancerBase::sendRecv(
     if(Pstream::parRun())
     {
 
-        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking, comm);
 
-        for(label i = 0; i < label(destinations.size()); ++i)
+        if (comm == UPstream::worldComm)
         {
-            UOPstream send(destinations[i], pBufs);
-            send << send_buffer[i];
+            for(label i = 0; i < label(destinations.size()); ++i)
+            {
+                UOPstream send(destinations[i], pBufs);
+                send << send_buffer[i];
+            }
+
+            pBufs.finishedSends();
+
+            ret.setSize(sources.size());
+            for(label i = 0; i < label(sources.size()); ++i)
+            {
+                UIPstream recv(sources[i], pBufs);
+                recv >> ret[i];
+            }
         }
-
-        pBufs.finishedSends();
-
-        ret.setSize(sources.size());
-        for(label i = 0; i < label(sources.size()); ++i)
+        else // now use Input inter-processor communications stream (UOPstream & UIPstream dose not work if use children communicator)
         {
-            UIPstream recv(sources[i], pBufs);
-            recv >> ret[i];
+            for(label i = 0; i < label(destinations.size()); ++i)
+            {
+                OPstream send(UPstream::commsTypes::scheduled, destinations[i], 0, 1, comm);
+                send << send_buffer[i];
+            }
+
+            ret.setSize(sources.size());
+            for(label i = 0; i < label(sources.size()); ++i)
+            {
+                IPstream recv(UPstream::commsTypes::scheduled, sources[i], 0, 1, comm);
+                recv >> ret[i];
+            }
         }
     }
 
diff --git a/src/dfCombustionModels/Make/options b/src/dfCombustionModels/Make/options
index fa21fbee..edeb1a3d 100644
--- a/src/dfCombustionModels/Make/options
+++ b/src/dfCombustionModels/Make/options
@@ -2,6 +2,7 @@
 
 EXE_INC = -std=c++14 \
     -Wno-old-style-cast \
+    $(PFLAGS) $(PINC) \
     $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \
     $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -17,12 +18,12 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/lagrangian/spray/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(DF_SRC)/dfCanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/functionObjects/field/lnInclude \
     -I$(CANTERA_ROOT)/include \
diff --git a/test/corrtest.cpp b/test/corrtest.cpp
index a2b68593..2099f3a4 100644
--- a/test/corrtest.cpp
+++ b/test/corrtest.cpp
@@ -25,8 +25,10 @@ float TGV400 = readTGV(1098,"2DTGV/4/data_T.xy");
 
 
 TEST(corrtest,df0DFoam_H2){
-    EXPECT_FLOAT_EQ(H2maxT,2588.48);   // compare the maximum temperature of H2 case 
-    EXPECT_FLOAT_EQ(H2midT,1021.41); // compare the temperature of H2 case at the maximum gradient when t = 0.000245s
+    // EXPECT_FLOAT_EQ(H2maxT,2588.48);   // compare the maximum temperature of H2 case 
+    // EXPECT_FLOAT_EQ(H2midT,1021.41); // compare the temperature of H2 case at the maximum gradient when t = 0.000245s
+    EXPECT_FLOAT_EQ(H2maxT,2586.21);   // compare the maximum temperature of H2 case 
+    EXPECT_FLOAT_EQ(H2midT,1020.71); // compare the temperature of H2 case at the maximum gradient when t = 0.000245s
 }