void user_source	(	size_t	dimension,
		size_t	num_points,
		double	y[],
		double	f[],
		double	growth[]
	)

Function that returns the growth and source terms.

Function where the user enters the evaluation of the source terms and growth terms.

External function to compute user defined source terms.

Author:: Rochan R. Upadhyay

Parameters:

dimension	: Dimension of the problem
num_points	: Number of quadrature points
y	: C style vector that contains quadrature points and weights. Details below.
f	: C style vector containing moments of the source terms of the Population Balance equation
growth	: C style vector that contains the growth terms (not the moments of the growth term)

The vector y contains the num_points weights W_i and the num_points * quadrature points x_{k,i}. With k = 0 to dimension-1 and i = 0 to n = 0 to num_points-1. The ordering is as follows: [ W_0,...,W_{num_points-1},x_{0,0},...,x_{0,num_points-1},x_{1,0},...,x_{1,num_points-1},...,x_{dimension-1,0},...,x_{dimension,num_points-1}]

For this example, the user should compute moments of the output variable and return them as source terms. The growth terms are unused and should be returned as zero. Note for GLD reconstruction only the first five moments and needed and hence only these moments are computed.

Date:: 10/12/2010

Computation of various quantities needed for the ASET model.

 double anum, adenom, term1, term2, time;
 double Z[num_points];
 anum=-0.21*(pow(g,0.5))*(pow(1.0-aLr,1.0/3.0));
 adenom=pow(rhoa*Ca*Ta*(pow(g,0.5)),1.0/3.0)*Ar;
 term1=anum/adenom;  
 term2=(1.0-aLc)/(rhoa*Ca*Ta*Ar);
 double Zf = Zf_factor*Z0;

Now loop over the quadrature points to compute the smoke height for each heat release rate. The vector y[] contains quadrature points and weights of the heat release rate distribution.

 for(size_t ipt=0; ipt < num_points; ipt++){
 Z[ipt] = Z0;

Evolve the ODE upto the critical time

 for( int it=0; it < num_steps; it++){
 time=(it+1)*dt;

Z[ipt]=Z[ipt]+dt*(term1*(pow(y[ipt+num_points],1.0/3.0)) (pow(Z[ipt]-Zf,5.0/3.0))-term2*y[ipt+num_points]); if(Z[ipt] <= Zf){Z[ipt]=Zf;}} }

Now compute the moments of Z at the critical time. Only the first five moments are required for the GLD reconstruction. The vector

 for (size_t l = 0; l < 6; l++){
     f[l] = 0.0;
     for (size_t i = 0; i < num_points; i++){
     f[l] += y[i]*pow(Z[i],(int)l);}}

Return control to the calling routines.

 return;    
 }