/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.


\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "fvOptions.H"
#include "simpleControl.H"
#include "OFstream.H"
#include <algorithm>


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Diffusion solver. "
    );
    #include "postProcess.H"
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
	#include "createFields_fi.H"

    
    Info<< "\nDiffusion equation\n" << endl;
    #include "geometrical_study.H"
	    scalar kk_eff_new=1; scalar kk_eff_old=1; scalar err_kk=10; // Definizione delle variabili
        scalarField fi_new; scalarField fi_old;  scalar FI_NEW=0; scalar FI_OLD=0;
		scalar err_fi=1; scalar err_fi_p2; scalar err_fi_p2_num=9; scalar err_fi_p2_den=0;
		scalar ii=0; scalar aa; 
		scalar PP; scalar PP_norm; scalar ee=3.204353268e-11;  // 200 MeV in Joule
        scalar volume=0;
		forAll(fi,cellI)
		{	
		volume += mesh_fi.V()[cellI];
		}
		 
        cout<<"Give a value for Power "; 	// Tramite questo e il comando sotto posso definire la otenxa da tastiera	
		cin>>PP; 
		
		
       // POWER METHOD
	   
		while (err_kk>1e-15 || err_fi>1e-13) // Il ciclo while del Power Method (PM) si trova sotto quello temporale in modo tale che OF salvi nella prima
		                                   // casella temporale ciò che ha prodotto questo ciclo while (OF salva le varibili in cartelle che crea ogni volta
										   // che si manda una simulazione)
        {
		fi_old=fi.internalField(); // copia del flusso attuale (quindi dell'iterazione precedente) in una variabile creata ad hoc
		kk_eff_old=kk_eff_new;	   // copia del k effettivo
        ii=ii+1;
		
        FI_OLD=0; // Azzeramento del numero dei neutroni prodotti (poiché uso un ciclo for sommerrei anche i neutroni delle iterazioni precedenti)
		FI_NEW=0;
 		
		forAll(fi,cellI) //calcolo del numero di neutroni prodotti nell'iterazione precedente
		{
		FI_OLD += (nisf[cellI]*fi[cellI]*mesh_fi.V()[cellI]);
		}
            
		fvScalarMatrix diff_term // Soluzione di L(fi)=Q dove Q è ipotizzato nel primo step 
		(   
	       
         -fvm::laplacian(DD,fi)
	     +fvm::Sp(sa,fi)  
		);
	    solve(diff_term == fvc::SuSp(one,qq)); 
				
				
		forAll(fi,cellI) // calcolo del numero di neutroni prodotti dopo L(fi)=Q quindi nell'iterazione attuale
		{
		FI_NEW += (nisf[cellI]*fi[cellI]*mesh_fi.V()[cellI]);
		}
    
		Info<<"Old neutrons "<<FI_OLD<<" new neutrons "<<FI_NEW<<endl; // stampa dei risultati
         

		
		// VALUTAZIONE DEL K EFFETTIVO
		
	    kk_eff_new=(kk_eff_old*FI_NEW)/FI_OLD; //calcolo di k effettivo
		
		err_kk=mag(kk_eff_new-kk_eff_old); // valutazione dell'errore su k effettivo
		 
		 
		// CALCOLO ERRORE SUL FLUSSO
		err_fi_p2_num=0;
		err_fi_p2_den=0; 
		 
		forAll(fi,cellI) // calcolo della norma - non è stata trovata una funziona adatta qundi è stato implementato un ciclo for. 
		{
		 err_fi_p2_num+=sqr(fi_old[cellI]-fi[cellI]); // numeratore
		 err_fi_p2_den+=(sqr(fi_old[cellI])); // denominatore
		}
		 		 
		err_fi=Foam::sqrt(err_fi_p2_num/err_fi_p2_den); // calcolo errore
		 
        Info<<"kk= "<<kk_eff_new<<" at iteration "<<ii<<" Error on kk "<<err_kk<<" Error on fi "<<err_fi<<nl<<endl;
		  
		  
		// CALCOLO DI Q CON IL FLUSSO DELL'ATTUALE ITERAZIONE	
		fvScalarMatrix qq_term
	    (
	     fvm::Sp(one,qq)        
        ); 
         
	     solve(qq_term == (1/kk_eff_new)*fvc::SuSp(nisf,fi));  
         // Q calcolato alla riga appena sopra è figlio del flusso dell'iterazione in cui ci si trova. Pertanto è Q(n+1) e quindi si usa kk_eff_new.
		 // L'operazione qui svolta è Q(n+1)=(1/kk(n))*F*fi(n). Se si usasse kk_eff_old l'equazione da risolvere sarebbe Q(n+1)=(1/kk(n-1))*F*fi(n)
        }
		// FINE DEL PM
		
		// NORMALIZZAZIONE DEL FLUSSO SECONDO LA POTENZA DATA IN INPUT
	    fi_new=fi.internalField(); // viene memorizzato il flusso proveniente dal PM in una variabile (il fatto che si chiami fi_new è solo perché era una 
		                           // varibile presente ma non ha alcun particolare significato)
		
	    forAll(fi,cellI) 
		{ 
         PP_norm+=nisf[cellI]*fi_new[cellI]*dx*dy*dz*ee; // calcolo della potenza prodotta
		}		 
		
		aa=PP/PP_norm; // calcolo del coefficiente da inserire davanti a fi per il calcolo del flusso affinché sia prodotta la potenza data in input
		
		forAll(fi,cellI) 
		{
		 fi[cellI]=aa*fi_new[cellI]; // ricalcolo del flusso con la costante aa. Ora il flusso se moltiplicato per la cross section di fissione, 
		                             // l'energia prodotta da fissione e integrato darà la potenza data in input
		}
		
		scalar stop;
		cout<<"I compute the flux. Do I stop? 1- Yes 2- No ";
		cin>>stop; 
		
		if (stop==1)
		{
		#include "operators.H"
		fi=fi/1e4;
		grad_fi_y=grad_fi_y/1e6;
		grad_DD_fi=grad_DD_fi/1e4;
		lap_fi=lap_fi/1e8;
		lap_dDD_fi=lap_dDD_fi/1e6;
		sa=sa/100;
		DD=DD*100;
		nisf=nisf/100;
		dDD=dDD*100;
		ds=ds/100;
		dnisf=dnisf/100;
		
		runTime++;
		
		fi.write();
		grad_fi_y.write(); 
		grad_DD_fi_y.write();
		lap_fi.write();
		lap_dDD_fi.write();
		DD.write();
		sa.write();
		nisf.write();
		dDD.write();
		ds.write();
		dnisf.write();
		
		exit (EXIT_FAILURE); 

		}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////GEOMETRICAL AND MESH PART/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



    #include "operators.H"
	scalar mash;
	cout<<"Please modify mesh: make autobuilding (1) or make your own (2) ";
    cin>>mash;

    #include "Make_mesh.H" 

   #include "createMesh.H"
   
	simpleControl simple(mesh);
    
	#include "createFields.H"

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////GEOMETRICAL AND MESH PART//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    Info<< "\nCalculating scalar transport\n" << endl;

    #include "CourantNo.H"
	
	scalar vv; 
    cout<<"Velocity = "; cin>>vv; vv=0.01*vv; 
	
	scalar delayed; 
	cout<<"Do you want consider delayed neutrons? 1- Yes 2- No ";
	cin>>delayed; 
if (delayed==1){	
#include "Weights_delayed.H"
	
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

            fvScalarMatrix AEqn
            (
               alfa*fvm::ddt(A)==
               gamma*fvm::div(phi,A)
              + delta*fvm::laplacian(A)
	          + fvm::Sp(eta,A)
			  + fvc::SuSp(lambda,gg)
             
            );
			
            AEqn.solve();
			A.write();
			
			aa=0;
			forAll(A,cellI) 
			{
			 aa=aa+A[cellI]*dy/HH;
			}
			forAll(A,cellI) 
			{
				AA[cellI]=aa; 
			}
			AA.write(); 
			
			
			fvScalarMatrix ggEqn{
			fvm::ddt(gg)==
			- fvm::Sp(lambda,gg)
			+ fvc::SuSp(zeta,A)
			};
            
			ggEqn.solve();
			gg.write(); 


		        aa=0;
		        forAll(A,cellI)
	                {  aa=aa+gg[cellI]*dy/HH;											                }
	               forAll(A,cellI)
		       {		
		        gg_int[cellI]=aa;
                        }
	                gg_int.write();
 			
    }

scalar coefficients; 
cout<<"Do I save coefficients? 1-Yes 2-No";
cin>>coefficients;
if(coefficients==1){
	alfa=0.01*alfa;
	delta=100*delta;
	eta=0.01*eta;
	zeta=0.01*zeta;
	
	alfa.write(); 
	gamma.write();
	delta.write();
	eta.write();
	zeta.write();
}
} else{
	
    
#include "Weights.H"
	
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        
        
            fvScalarMatrix AEqn
            (
               alfa*fvm::ddt(A)==
               gamma*fvm::div(phi,A)
              + delta*fvm::laplacian(A)
	          + fvm::Sp(eta,A)
             
            );
			
            AEqn.solve();
			A.write();
			
			aa=0;
			forAll(A,cellI) 
			{
			 aa=aa+A[cellI]*dy/HH;
			}
			forAll(A,cellI) 
			{
				AA[cellI]=aa; 
			}
			AA.write(); 
              
    }

scalar coefficients; 
cout<<"Do I save coefficients? 1-Yes 2-No ";
cin>>coefficients;
if(coefficients==1){
	alfa=0.01*alfa;
	delta=100*delta;
	eta=0.01*eta;
	
	alfa.write(); 
	gamma.write();
	delta.write();
	eta.write();
}

}


    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //
