/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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
    (
        "Passive scalar transport equation solver." // Da ignorare
    );
	
	
    #include "postProcess.H"
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    //#include "createMesh.H"
    #include "createFields_fi.H"
	
    //simpleControl simple(mesh_fi);

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

    Info<< "\nDiffusion equation for two groups\n" << endl;
// N.B.
// Il metodo delle potenze nel caso multigruppo prevede come semrpe L(fi)=Q e F(fi)=Q dove L e F diventano delle matrici e fi un vettore fatto dai diversi 
// flussi energetici. Qui non si è riuscito ad implementare le matrici L e F e si è fatta una divisione "non convenzionale" in L1 e L2, Fi e F2, Q1 e Q2.
// Essi corrispondono ai diversi flussi energetici (1 ad alta energia, 2 a bassa energia perciò) e al termine di diffusione, fissione e sorgente. Per esempio
// L2 è la parte di diffusione del flusso di neutroni a bassa energia e corrisponde a -lap(DD2,fi2) + Sigma_abs_2*fi2 - Sigma_1->2*fi1. 
// Si noti che: in L1 compare solo fi1; in F1 compare anche il termine fi2 e di ciò se ne tiene conto andando a risolvere prima l'equazione per fi1
// poi per fi2 e poi si usano questi per calcolare Q1; F2 è 0 e pertanto il termine Q2 è sempre 0 e non necessita di ricalcolo
 
	    // DICHIARAZIONE DELLE VARIABILI
		#include "Geometrical_study.H"
    	scalar mmm=0; scalar mm=0;scalar ii=0;
		scalar kk_eff_old=1; scalar kk_eff_new=1;scalar err_kk=10;
		scalar FI_NEW=0; scalar FI_OLD=0; scalarField fi1_new; scalarField fi1_old; scalarField fi2_new; scalarField fi2_old;
        scalar err_fi;scalar err_fi_p2_num=9; scalar err_fi_p2_den;
		scalar PP_norm; scalar PP=1000;cout<<"\nGive a value for Power  "; cin>>PP;// Questo e il comando successivo permettono di poter prendere la potenza come input da tastiera	
		scalar volume=0;scalar aa; scalar ee=3.204353268e-11; scalar nn_cell=0; // 200 MeV in Joule      
			
		forAll(fi1,cellI) // Calcolo del volume
		{	
		volume += mesh_fi.V()[cellI];
		nn_cell=cellI;
		}
		// VARIABILI PER L'IMPORTANZA
        scalar jj=0; scalar bb; 
		scalar kk_eff_psi_old=1;scalar kk_eff_psi_new=1;scalar err_kk_psi=10;
		scalar PSI_NEW=0;scalar PSI_OLD=0; scalarField psi1_new; scalarField psi1_old; scalarField psi2_new; scalarField psi2_old; 
		scalar err_psi; scalar err_psi_p2_num=9; scalar err_psi_p2_den;
		scalar PP_norm_psi; scalar PP_psi=1; cout<<"Give a value for normalization "; cin >> PP_psi;
	    // CONVERSIONE DA n/cm^2/s a n/m^2/s 
		
        //POWER METHOD
		while (err_kk>1e-15 || err_fi>1e-13) // Il ciclo while del Power Method (PM) si trova dentro 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, riconducibili al timestep usato )
        {
		fi1_old=fi1.internalField(); // memorizzazione dei vecchi flussi ossia quelli dell'iterazione precedente
        fi2_old=fi2.internalField();		
			
		kk_eff_old=kk_eff_new;	// memorizzazione del precedente k effettivo
        ii=ii+1;				    
        FI_OLD=0; // annullamento delle variabii per la memorizzazione del numero di neutroni prodotti
		FI_NEW=0;
		
		forAll(fi1,cellI) // calcolo del numero di neutroni prodotti nell'iterazione precedente. Sarebbe un integrale nello spazio e nell'energia ma non essendoci
		                  // una vera e propria divisione dell'energia si somma senza alcun "dE"
		{
		FI_OLD += (nisf1[cellI]*fi1[cellI]*mesh_fi.V()[cellI]+nisf2[cellI]*fi2[cellI]*mesh_fi.V()[cellI]);
		}
         
        // SOLUZIONE EQUAZIONI
		
		fvScalarMatrix diff_term_1  // soluzione di L1(fi_1) = Q1 Si noti che compare solo il termine fi1 in questa parte di equazione. Q1 che nel caso del         
		                            // primo step è ipotizzato, ha in sé sia fi_1 che fi_2. L'operazione qui presente è L1(fi_1(n)) = Q1(n)
		(   
	       
         -fvm::laplacian(DD1,fi1)
	     +fvm::Sp(sa1,fi1)  
		);

	    solve(diff_term_1 == fvc::SuSp(one,qq1)); 
		
		fvScalarMatrix diff_term_2 // soluzione di L2(fi_2)=Q2 Si noti che questa volta vi è anche il termnie fi1 che è presente nell'equazione ma come termine 
		                           // noto. Questo è quello appena calcolato con l'equazione sovrastante. Si risolve dunque L(fi_2(n))=Q2(n) dove in L2 (ma poi 
								   // OF lo risolve come se fosse in Q2?) vi è anche fi_1(n)
								   
		(   
	      
         -fvm::laplacian(DD2,fi2)
	     +fvm::Sp(sa2,fi2)
		 -fvc::SuSp(ss12,fi1)		 
		);

	    solve(diff_term_2 == fvc::SuSp(one,qq2)); // Q2 che qui si vede è posto uguale a 0 sempre
	
	
        // CALCOLO DEL NUMERO DI NEUTRONI PRODOTTI CON IL NUOVO FLUSSO
		forAll(fi1,cellI) // calcolo dei neutroni prodotti con i nuovi flussi
		{
		FI_NEW += (nisf1[cellI]*fi1[cellI]*mesh_fi.V()[cellI]+nisf2[cellI]*fi2[cellI]*mesh_fi.V()[cellI]);
		}
    
            
	    kk_eff_new=(kk_eff_old*FI_NEW)/FI_OLD; // calcolo del nuovo k effettivo
		
			
			
		// SOLUZIONE DELL'EQUAZIONE PER Q1	
		fvScalarMatrix qq_term_1 // calcolo di Q1 come F(fi(n))=Q1 si noti che fi(n) contiene sia fi_1 che fi_2. Scrivendo le equazioni per la diffusione
		                         // a 2 gruppi si trovano L e F e che F è una matrice che ha una riga fatta dalle due cross section di fissione e una riga                
								 // fatta di 0 (per il secondo gruppo). Ne viene che L(fi_2)=0 e che pertanto non occorre nessun ricalcolo di Q2
	    (
	     fvm::Sp(one,qq1)        
        ); 
      
	     solve(qq_term_1 == (1/kk_eff_new)*(fvc::SuSp(nisf1,fi1)+fvc::SuSp(nisf2,fi2)));
		 		
		
		 err_kk=mag(kk_eff_new-kk_eff_old); // calcolo dell'errore su k
		 
		 
		// CALCOLO DELL'ERRORE SU FI - CALCOLO DELLA NORMA 
		err_fi_p2_num=0; //azzeramento del numeratore per il calcolo della norma
		err_fi_p2_den=0; //azzeramento del denominatore per il calcolo della norma
		 
		forAll(fi1,cellI) // calcolo della norma - non è stata trovata una funziona adatta qundi è stato implementato un ciclo for. 
		{
		 err_fi_p2_num+=sqr(fi1_old[cellI]-fi1[cellI]); 
		 err_fi_p2_num+=sqr(fi2_old[cellI]-fi2[cellI]);
		 err_fi_p2_den+=(sqr(fi1_old[cellI]));
		 err_fi_p2_den+=(sqr(fi2_old[cellI]));
		}
		 		 
		err_fi=Foam::sqrt(err_fi_p2_num/err_fi_p2_den); 
         
		 
		Info<<"kk= "<<kk_eff_new<<" Error on kk "<<err_kk<<" Error on fi "<<err_fi<<" at "<<ii<<endl; // stampa dei risultati
        }
		// FINE DEL CICLO WHILE PER IL FLUSSO 
		
		fi1_new=fi1.internalField(); // memorizzazione del flusso calcolato col PM in una variabile temporale
		fi2_new=fi2.internalField();
		
	    forAll(fi1,cellI) // calcolo della potenza prodotta
		{ 
         PP_norm+=(nisf1[cellI]*fi1_new[cellI]*dx*dy*dz*ee+nisf2[cellI]*fi2_new[cellI]*dx*dy*dz*ee);
		}		 
		
		aa=PP/PP_norm; // calcolo della costante aa per la normalizzazione del flusso secondo la potenza scelta
		
		forAll(fi1,cellI) // aggiornamento del flusso con la nuova costante davanti. 
		{
		 fi1[cellI]=aa*fi1_new[cellI]; 
		 fi2[cellI]=aa*fi2_new[cellI];
		}
		
		Info<<"\n\n\nkk= "<<kk_eff_new<<" Error on kk "<<err_kk<<" Error on fi "<<err_fi<<endl; // stampa dei risultati

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////		
if(mm==0)
{
mm=mm+1;	
cout<<"\n\n\n\n//////////////////////////////////////////      ADJOINT    //////////////////////////////////////////\n\n\n\n";        
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////		
		// INIZIO CICLO WHILE PER L'AGGIUNTO
		
		while (err_kk_psi>1e-15 || err_psi>1e-13)
        {
		kk_eff_psi_old=kk_eff_psi_new;	//memorizzazione del vecchio k effettivo, dell'iterazione precedente
        jj=jj+1;
		psi1_old=psi1.internalField(); //memorizzazione delle vecchie importanze
        psi2_old=psi2.internalField();		
        PSI_OLD=0;
		PSI_NEW=0;
		
		forAll(psi1,cellI) // calcolo dell'importanza totale generata con il valore precedente
		{
		PSI_OLD += (nisf1[cellI]*psi1[cellI]*mesh_fi.V()[cellI]+nisf2[cellI]*psi1[cellI]*mesh_fi.V()[cellI]);
		}
         
		//SOLUZIONE DELLE EQUAZIONI 
		 fvScalarMatrix diff_term_psi_2 //soluzione di L2(fi2)=Q2 dove in Q2 è nascosto il termine fissione applicato a fi1. Nel caso di primo step è ipotizzato
		(   
	      
         -fvm::laplacian(DD2,psi2)
	     +fvm::Sp(sa2,psi2)
		  
		);

	    solve(diff_term_psi_2 == fvc::SuSp(one,qq_psi_2));
   		
		fvScalarMatrix diff_term_psi_1 //soluzione di L1(fi1)=Q1 dove Q1 è il termine di sorgente. Si noti che questa volta compare fi2 che dovrebbe essere 
                                   // considerato come Q1 come termine noto. Il termine fvc::SuSp(ss12,fi2) fa riferimento al passaggio di neutroni 
                                   // da un gruppo all'altro e fi2 è quello appena calcolato nell'equazione sovrastante								   
		(   
	      
         -fvm::laplacian(DD1,psi1)
	     +fvm::Sp(sa1,psi1)  
		 - fvc::SuSp(ss12,psi2)
		);

	    solve(diff_term_psi_1 == fvc::SuSp(one,qq_psi_1)); 
		
		 
		forAll(psi1,cellI) // calcolo dell'importanza totale prodotta con le due nuove importanza
		{
		PSI_NEW += (nisf1[cellI]*psi1[cellI]*mesh_fi.V()[cellI]+nisf2[cellI]*psi1[cellI]*mesh_fi.V()[cellI]);
		}
    
            
	    kk_eff_psi_new=(kk_eff_psi_old*PSI_NEW)/PSI_OLD; // calcolo di k effettivo
		
		
		//CALCOLO DI Q1 E Q2
		
		fvScalarMatrix qq_term_psi_1 //Per Q1 si usa fi1 appena trovato quindi l'operazione dovrebbe essere Q1(n+1)=F1(fi1(n))
	    (
	     fvm::Sp(one,qq_psi_1)        
        ); 
      
	     solve(qq_term_psi_1 == (1/kk_eff_psi_new)*(fvc::SuSp(nisf1,psi1)));
		 		
		fvScalarMatrix qq_term_psi_2 // Questa volta il termine Q2 non è più zero ma deve essere riaggiornato con il flusso fi1 prima calcolato
	    (
	     fvm::Sp(one,qq_psi_2)        
        ); 
      
	     solve(qq_term_psi_2 == (1/kk_eff_psi_new)*fvc::SuSp(nisf2,psi1)); 
		
        // VALUTAZIONE DEGLI ERRORI SU K E FI		
		err_kk_psi=mag(kk_eff_psi_new-kk_eff_psi_old); // calcolo errore su k
		 
		err_psi_p2_num=0; // azzeramento del numeratore per il calcolo della norma quindi per l'errore su fi
		err_psi_p2_den=0; // azzeramento del denominatore per il calcolo della norma quindi per l'errore su fi
		 
		forAll(psi1,cellI) // calcolo della norma - non è stata trovata una funziona adatta qundi è stato implementato un ciclo for. 
		{
		 err_psi_p2_num+=sqr(psi1_old[cellI]-psi1[cellI]); 
		 err_psi_p2_num+=sqr(psi2_old[cellI]-psi2[cellI]); // Si noti che le due sommatorie, al numeratore e al denominatore, hanno come addendi 
                                                        // sia le differenze di fi1 sia quelle di fi2 (es):
														// (fi_1_old[i]-fi_1_new[i])^2 + ... + (fi_2_old[i]-fi_2_new[i])		 
		 err_psi_p2_den+=(sqr(psi1_old[cellI])); 
		 err_psi_p2_den+=(sqr(psi2_old[cellI]));
		}
		 		 
		err_psi=Foam::sqrt(err_psi_p2_num/err_psi_p2_den); 
         
		 
		 Info<<"kk= "<<kk_eff_psi_new<<" Error on kk "<<err_kk_psi<<" Error on psi "<<err_psi<<" at "<<ii<<endl; //Stampa dei risultati 
        }
		
		// NORMALIZZAZIONE DELL'AGGIUNTO SECONDO LA POTENZA DATA IN INPUT
        	
		fi1_new=fi1.internalField(); 
		fi2_new=fi2.internalField();
		psi1_new=psi1.internalField(); 
		psi2_new=psi2.internalField();
		
	    forAll(psi1,cellI) // calcolo della potenza prodotta
		{ 
         PP_norm_psi+=(fi1_new[cellI]*psi1_new[cellI]*dx*dy*dz+fi2_new[cellI]*psi2_new[cellI]*dx*dy*dz);
		}		 
		
		bb=PP_psi/PP_norm_psi; // calcolo della costante aa per la normalizzazione del flusso secondo la potenza scelta
		
		forAll(psi1,cellI) // aggiornamento del flusso con la nuova costante davanti. 
		{
		 psi1[cellI]=bb*psi1_new[cellI]; 
		 psi2[cellI]=bb*psi2_new[cellI];
		}
		
		Info<<"\n\n\nkk= "<<kk_eff_psi_new<<" Error on kk "<<err_kk_psi<<" Error on fi "<<err_psi<<endl; //Stampa dei risultati 

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////SALVATAGGIO DATI ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
scalar zero=0; 

scalar fermati;
cout<<"I computed flux. Do I stop? 1-Yes 2-No ";
cin>>fermati; 
if (fermati==1)
    { 	

    #include "operators.H"
    
   	fi1=fi1/1e4;
	fi2=fi2/1e4;
	psi1=psi1/1e4;
	psi2=psi2/1e4;
	grad_DD1_fi1_y=grad_DD1_fi1_y/1e4;
	grad_DD2_fi2_y=grad_DD2_fi2_y/1e4;
	lap_fi1=lap_fi1/1e8;
    lap_fi2=lap_fi2/1e8;
	lap_dDD1_fi1=lap_dDD1_fi1/1e6;
    lap_dDD2_fi2=lap_dDD2_fi2/1e6;
    DD1=DD1*100;
    DD2=DD2*100;
	dDD1=dDD1*100;
    dDD2=dDD2*100;
    sa1=sa1/100;
    sa2=sa2/100;
    ds1=ds1/100;
    ds2=ds2/100;
    nisf1=nisf1/100;
    nisf2=nisf2/100;
    ss12=ss12/100;
    ds12=ds12/100;

    runTime++;
	
	DD1.write();
	DD2.write();
	dDD1.write();
	dDD2.write();
    sa1.write();
    sa2.write();
    ds1.write();
    ds2.write();
    nisf1.write();
    nisf2.write();
    dnisf1.write();
    dnisf2.write();
    ds12.write();
     fi1.write(); 
	 fi2.write(); 
	 psi1.write(); 
	 psi2.write();
     ss12.write();
     ds12.write();
	 grad_fi1_y.write(); 
	 grad_fi2_y.write(); 
	 grad_DD1_fi1_y.write();
	 grad_DD2_fi2_y.write(); 
	 lap_fi1.write();
	 lap_fi2.write();
	 lap_dDD1_fi1.write();
	 lap_dDD2_fi2.write();


	
	 exit (EXIT_FAILURE); 
	}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////// CALCOLO DEI PESI ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

cout<<"\n\n\n\n//////////////////////////////////////////    WEIGHT COMPUTATION   //////////////////////////////////////////\n\n\n\n";        
cout<<"Modify mesh to compute the amplitude. 1 for automatic process 2 if you prefer to do yourself ";
scalar mash;    
cin>>mash;	

#include "Make_mesh.H"  
	
// Inizialization of weights
scalar sum; scalar kk;
#include "operators.H"
#include "createMesh.H"


simpleControl simple(mesh);

#include "createFields.H"


       scalar vv1; scalar vv2;
       cout<<"Velocity 1= "; cin>>vv1;  cout<<"Velocity 2= ";cin>>vv2;
	   vv1=0.01*vv1; vv2=0.01*vv2;


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////// CALCOLO AMPIEZZE ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


cout<<"\n\n\n\n//////////////////////////////////////////    Amplitude computation   //////////////////////////////////////////\n\n\n\n";   
scalar delayed;	
cout<<"Delayed? 1-Yes 2-No "; 
cin>>delayed;
scalar salvare_dati; 
cout<<"Save weights? 1-Yes 2-No "; 
cin>>salvare_dati;


if (delayed==2)
{	
    #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();     
              
			sum=0; 
			forAll(A,cellI)
            { 
			 sum=sum+A[cellI]*dy/HH;
	        }
			 A_int=sum;
             A_int.write();


		runTime.write();
    }
if (salvare_dati==1)
		{
	   alfa=alfa/100; 	
       delta=delta*100;	   
	   eta=(eta/100);

	   alfa.write();	   	
	   gamma.write();
	   delta.write();
	   eta.write();
	   
		}

}
else
{
#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();
            
			sum=0; 
			forAll(A,cellI)
            { 
			 sum=sum+A[cellI]*dy/HH;
	        }
			 A_int=sum;
			
			A.write(); 
			A_int.write();

			  					
		
		fvScalarMatrix ggEqn
            (
               fvm::ddt(gg)==
			   -fvm::Sp(lambda,gg)
			   +fvc::SuSp(zeta,A)

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

		  forAll(A,cellI)
	           {
		  gg_int[cellI]=aa;
		  }
	        gg_int=aa;
	        gg_int.write();


		runTime.write();
		

	
   }


	   if (salvare_dati==1)
		{
	   alfa=alfa/100; 	   
	   delta=(delta*100);
	   eta=(eta/100);
	   zeta=(zeta/100);
	   
	   alfa.write();	   
	   gamma.write();
	   delta.write();
	   eta.write();
	   zeta.write();
	   
		}




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

    return 0;
	}


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