From d0bf876a7cefdf14f1d39f2ea69ade7274d41f89 Mon Sep 17 00:00:00 2001
From: michael <michael@localhost.localdomain>
Date: Fri, 10 Dec 2021 14:34:48 +0100
Subject: [PATCH] Tested mesh fast equivalence when giving an input field + MPI
 parallelisation of the class

---
 CoreFlows/Models/inc/DiffusionEquation.hxx | 14 ++-----
 CoreFlows/Models/src/DiffusionEquation.cxx | 44 ++++++++++++++++++++--
 2 files changed, 44 insertions(+), 14 deletions(-)

diff --git a/CoreFlows/Models/inc/DiffusionEquation.hxx b/CoreFlows/Models/inc/DiffusionEquation.hxx
index ec2692c..16e24c5 100755
--- a/CoreFlows/Models/inc/DiffusionEquation.hxx
+++ b/CoreFlows/Models/inc/DiffusionEquation.hxx
@@ -106,17 +106,9 @@ public :
 	vector<string> getInputFieldsNames();
 	void setInputField(const string& nameField, Field& inputField );//supply of a required input field
 
-	void setFluidTemperatureField(Field coupledTemperatureField){
-		coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
-		_fluidTemperatureField=coupledTemperatureField;
-		_fluidTemperatureFieldSet=true;
-	};
-	void setFluidTemperature(double fluidTemperature){
-	_fluidTemperature=fluidTemperature;
-	}
-	Field& getFluidTemperatureField(){
-		return _fluidTemperatureField;
-	}
+	void setFluidTemperatureField(Field coupledTemperatureField);
+	void setFluidTemperature(double fluidTemperature){	_fluidTemperature=fluidTemperature;	}
+	Field& getFluidTemperatureField(){  return _fluidTemperatureField;	}
 	
 	/*** get output fields names for postprocessing or coupling ***/
 	vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
diff --git a/CoreFlows/Models/src/DiffusionEquation.cxx b/CoreFlows/Models/src/DiffusionEquation.cxx
index 4c4a718..f769870 100755
--- a/CoreFlows/Models/src/DiffusionEquation.cxx
+++ b/CoreFlows/Models/src/DiffusionEquation.cxx
@@ -269,10 +269,35 @@ void DiffusionEquation::initialize()
 	//Linear solver
 	KSPCreate(PETSC_COMM_WORLD, &_ksp);
 	KSPSetType(_ksp, _ksptype);
-	// if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
-	KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 	KSPGetPC(_ksp, &_pc);
-	PCSetType(_pc, _pctype);
+	if(_mpi_size==1 )
+		PCSetType(_pc, _pctype);
+	else
+	{
+		PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
+		if(_pctype != (char*)&PCILU)//Default pc type is ilu
+		{
+			PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
+			PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
+			//If the above setvalue does not work, try the following
+			/*
+			KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
+			KSP * subKSP;
+			PC subpc;
+			int nlocal;//nb local blocs (should equal 1)
+			PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
+			if(nlocal==1)
+			{
+				KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
+				KSPGetPC(subKSP[0],&subpc);
+				PCSetType(subpc,_pctype);
+			}
+			else
+				throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
+			*/ 
+		}
+	}
+	KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 
 	_initializedMemory=true;
 	save();//save initial data
@@ -890,6 +915,9 @@ DiffusionEquation::getOutputField(const string& nameField )
 void
 DiffusionEquation::setInputField(const string& nameField, Field& inputField )
 {
+	if(!_initialDataSet)
+		throw CdmathException("!!!!!!!! DiffusionEquation::setInputField set initial field first");
+
 	if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
 		return setFluidTemperatureField( inputField) ;
 	else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
@@ -900,3 +928,13 @@ DiffusionEquation::setInputField(const string& nameField, Field& inputField )
         throw CdmathException("DiffusionEquation::setInputField error : Unknown Field name");
     }
 }
+
+void 
+DiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
+	if(!_initialDataSet)
+		throw CdmathException("!!!!!!!! DiffusionEquation::setFluidTemperatureField set initial field first");
+
+	coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
+	_fluidTemperatureField=coupledTemperatureField;
+	_fluidTemperatureFieldSet=true;
+};
-- 
2.39.2