/** 
# Intro 

This is a minimal example for implementing a 'fan' after geometry examples that can be found [here](http://basilisk.fr/sandbox/Antoonvh/puck.c). The mentioned function is updated to be applicable for general phi and theta. Also this script adds a velocity based on a stated fan power. All functions can be found at the end of this file.

# Power to velocity forcing and expected exit velocity
The implementation is based on a momentum forcing. An amount of work per time is done on the flow with a certain power, P. Using the work-energy theorem, which states that the work done by all acting forces on a particle is equal the change in kinetic energy, this can be expressed as $\Delta E_{kin}=P\Delta t$. This directly implies that,

$$\vec{v}^2_{t+\Delta t} - \vec{v}^2_{t} = \frac{2\textrm{P}\Delta t}{m}\hat{n}^2$$

which in turn implies for every spatial dimension,

$$ v^2_{i,t+\Delta t} - v^2_{i,t} = \frac{2\textrm{P}\Delta t}{m} n_i^2 $$

From these equations the direction of $\hat{n}$ is not yet defined and $|\hat{n}|=1$. When it is assumed that the force applied to the flow is normal to the plane in which the fan blades rotate then $\hat{n}$ is the normal vector of this plane. Rearranging and taking the square root results in,
$$	\vec{v}_{i,t+\Delta t} = \sqrt{\vec{v}^2_{i,t} + \frac{2\textrm{P}\Delta t}{m}\hat{n}^2_i } $$

To estimate the velocity to which air is accelerated by this forcing we assume that it is applied over a finite time from $t_{start}$ to $t_{end}$, moving a distance $w$ in the process:
$$		w=\int_{t_{start}}^{t_{end}} v\left(t\right)\textrm{dt}$$
Assuming a zero entrance velocity, $v$ is given at every instance as:
$$\textrm{E}_{\textrm{kin}}=\textrm{P}\Delta t=0.5mv^2$$
$$v\left(t\right)=\sqrt{2\textrm{P}t/m}$$
P, $w$ and m are set by the user thus this leaves two equations and two unknowns. Solving for the exit velocity gives:
$$	v_{end} = \left(\frac{3\textrm{P}w}{m}\right)^{1/3}$$


# The code
First we include relevant files for 2D navier stokes solving and geometric fun.
*/
//#include "grid/octree.h" 		// For 3D --> to demanding for server
#include "navier-stokes/centered.h"     // Navier stokes 
#include "fractions.h"	      		// Needed to compute fan fractions

/** Required functions and structures */
void init_rotor();
void rotor_update();
void rotor_coord();
void rotor_forcing();

struct sRotor {	
	double rampT;		// Time to start up rotor
	double P, Prho;		// Power, powerdensity 
	double R, W, A, V;	// Diameter, Thickness, Area ,Volume
	double diaVol;		// Diagnosed rotor volume 
	double x0, y0, z0;	// Origin of rotor
    	double xt, yt, zt;      // Translation of origin per rotate event
        double theta, phi;	// Polar and Azimuthal angle 
   	double thetat, phit;    // Addition to angles per rotate event
   	double Work;            // Work done by rotor
   	double cu;              // Characteristic velocity
   	bool rotate;            // Rotation yes or no
	coord nf, nr;		// Normal vector fan, rotation 
};

/** Some global vars */
int minlevel, maxlevel;         // Grid depths
double eps;			// Adaptivity criterion
struct sRotor rot;  		// Rotor details structure 
scalar fan[];			// Fan volume fraction
scalar lev[]; 			// for movie diagnostic level 
scalar Ekin[];  		// for movie kinetic energy 

/** Starting with the most loved c function */

int main() {	
	minlevel = 3;
  	maxlevel = 8;

   	L0 = 10.;
   	X0 = Y0 = Z0 = 0.;

    	init_grid(1<<7);

	foreach_dimension(){
		periodic (left);
		u.x.refine = refine_linear; 		// Momentum conserved
        }

	fan.prolongation = fraction_refine;		// Fan is a volume fraction
	p.refine = p.prolongation = refine_linear;

	DT = 10E-5;					// For poisson solver 
        TOLERANCE=10E-6;				// For poisson solver 
	CFL = 0.5;					// CFL condition

    	run();						// Start simulation 	
}

/** Initialisation */
event init(t = 0){
	rot.rotate = true; 				// Turn on rotation
	init_rotor();					// See init function for details 
	fan.prolongation = fraction_refine;		// Tell basilisk that fan is volumetric field
	refine (fan[] > 0. && level < maxlevel); 	// Refine where fan is
	eps = 0.07*rot.cu;				// Wavelet estimator based on flow velocity
}

/** Reset to defaults for poisson solver */
event reset_pois(i = 10){
  DT = 1;
  TOLERANCE = 10E-3;
}

/** Adaptivity */
event adapt(i++) {
	adapt_wavelet((scalar *){fan,u},(double []){0.,eps,eps},maxlevel,minlevel);
}

/** Forcing by the rotor */
event forcing(i = 1; i++) {
	rotor_coord();
	rotor_forcing();	
}

/** Rotate the rotor */
event rotate(t+=0.1) {
    if(rot.rotate) { 
        // Change center  
        rot.x0 += rot.xt;
        rot.y0 += rot.yt;
        rot.z0 += rot.zt;

        // Change angles 
        rot.theta += rot.thetat;
        rot.phi += rot.phit;

        rotor_update();
    }
}

/** Visual output */
event movies(t+=1){
    fprintf(stderr, "t=%g\n", t);
    foreach(){
   	lev[] = level;
 	Ekin[] = sqrt(sq(u.y[]) + sq(u.x[]));
    }
    boundary({lev, Ekin});

    output_ppm (lev, file = "ppm2mp4 level.mp4", n = 512, linear = false, min = minlevel, max = maxlevel);
    output_ppm (Ekin, file = "ppm2mp4 ekin.mp4", n = 512, linear = true);


}

/** End the simulation */
event end(t=180){
}
/** 
# Results
![Kinetic energy density](fan/ekin.mp4)
![Grid depth](fan/level.mp4)
*/

/** 
# Functions 
*/

/** Function returning the sRotor structure, includign default properties */
void init_rotor() {
	rot.Work = 0.;
    	if(!rot.rampT)
	    	rot.rampT = 1.; // fan at full thrust after 1 second, linear ramp
    	if(!rot.R)
	    	rot.R = L0/80.;     
    	if(!rot.W)
	    	rot.W = rot.R/5;    
    	if(!rot.Prho)                  
     		rot.Prho = L0/10.;		
   	if(!rot.x0)
    		rot.x0 = L0/2.;
    	if(!rot.y0)
	    	rot.y0 = L0/2.;
    	if(!rot.z0){
        #if dimension == 2
            	rot.z0 = 0.;
        #elif dimension == 3
            	rot.z0 = L0/2.;
        #endif
        }
    	if(!rot.theta)
	    	rot.theta = 90.*M_PI/180.;      // Polar angle
    	if(!rot.phi)
	    	rot.phi = 90.*M_PI/180.;        // Azimuthal angle 
  
    	if(rot.rotate) {
                // determines for every rotate event how much needs to be added to x, y, z, tehta, phi.
        	rot.xt = 0;
        	rot.yt = 0;
        	rot.zt = 0;
        	rot.thetat = 0.;
       	 	rot.phit = -0.1*M_PI/180.; // note that this rate is linked to the frequency of the rotate event
    	} else {
       		rot.xt = 0;
        	rot.yt = 0;
        	rot.zt = 0;
        	rot.thetat = 0.;
        	rot.phit = 0.;
    }
	rotor_update();
}

/** Updating relevant rotor variables */
void rotor_update() {
   	// Normal vectors to the fan plane, note that z and y are turned around for consistency regarding polar coordinates (basilisk uses y for the vertical)
    	rot.nf.x = sin(rot.theta)*cos(rot.phi);
	rot.nf.y = sin(rot.theta)*sin(rot.phi);
	rot.nf.z = cos(rot.theta);

	rot.nr.x = sin(rot.theta)*cos(rot.phi);
   	rot.nr.y = sin(rot.theta)*sin(rot.phi);
    	rot.nr.z = cos(rot.theta);

   	#if dimension == 2	
		rot.A = 2*rot.R*rot.W;
	#elif dimension == 3    	
		rot.A = sq(rot.R)*M_PI;      
	#endif

        // Volume is slice of a sphere
  	#if dimension == 2
		rot.V = 1.*rot.A;
	#elif dimension == 3 
		rot.V = 4.*M_PI*pow(rot.R,3.)/3. - 
			2*M_PI*pow(rot.R-rot.W/2., 2.)/3.*(2*rot.R + rot.W/2.);
	#endif

	rot.P = rot.V*rot.Prho;
    	rot.cu = pow(3*rot.Prho*rot.W, 1./3.);
}

/** Function returning the volume fractions of a fan object */
void rotor_coord() {
    	scalar sph[], plnu[], plnd[];
   	fraction(sph, -sq((x - rot.x0)) - sq((y - rot.y0)) - sq((z - rot.z0)) + sq(rot.R)); // sphere
  	fraction(plnu,  rot.nr.x*(x - rot.x0) + rot.nr.y*(y - rot.y0) + rot.nr.z*(z - rot.z0) + rot.W/2.); // upper plane
   	fraction(plnd, -rot.nr.x*(x - rot.x0) - rot.nr.y*(y - rot.y0) - rot.nr.z*(z - rot.z0) + rot.W/2.); // lower plane

	foreach () {
    		fan[] = sph[]*plnu[]*plnd[]; // estimate for the overlappnig volumes, not that this is not exact at the boundaries of the 'fan', this is corrected for later in the focring by comparing diagnosed to theoretical volume. 
   	}
	boundary({fan}); // update boundaries
}

/** Function returning new velocities based on a rotor forcing.
This is a function of powerdensity, width, direction, ramp-time, and diagnosed volume */
void rotor_forcing(){
	double tempW = 0.;
	double w, wsgn, damp, usgn, utemp, corrP;
	foreach(reduction(+:tempW)) {		
		if(fan[] > 0.) {
			foreach_dimension() {
			wsgn = sign(rot.nf.x*u.x[]) + (sign(rot.nf.x*u.x[]) == 0)*sign(rot.nf.x); // Get the direction of the flow relative to forcing direction 
			damp = rot.rampT > t ? t/rot.rampT : 1.;                                  // Linear ramp for starting up the forcing
			corrP = rot.diaVol > 0. ? rot.V/rot.diaVol : 1.;                         // Correction if diagnosed volume differs from theoretical volume
			w = wsgn*fan[]*damp*sq(rot.nf.x)*(2./rho[])*(corrP*rot.P/rot.V)*dt;      // additional kinetic energy
			tempW += 0.5*rho[]*w*dv(); // save work done to asses total work in diagnostics

			// New kinetic energy
			utemp = sq(u.x[]) + w;
                        // New sign of the velocity
			usgn = 	  1.*(u.x[] >= 0)*(utemp > 0) +
			    	 -1.*(u.x[] >= 0)*(utemp < 0) +
		 		  1.*(u.x[] <  0)*(utemp < 0) +
				 -1.*(u.x[] <  0)*(utemp > 0); 
                        // New velocity
			u.x[] = usgn*sqrt(fabs(utemp)); 
		}
		}
	}
	rot.Work += tempW; // Diagnostics
}