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}]

The user should compute moments of the source terms aand return the vector f of size num_points * (dimension + 1). The growth terms which are velocities in internal space, should be provided in the vector growth again of size num_points * (dimension + 1).

Date:: 10/12/2010

Definition of the size of the problem itot, i.e. the number of moments that need to be evolved. In this example, dimension = 1

 size_t itot = dimension*num_points + num_points;

Outer loop over the total number of moments.

 for( size_t i = 0; i < itot; i++ )
 {

Computing source terms for the moment equations. For QMOM the source terms involve integrals of the term that models growth and diffusion in internal space. The growth terms are given by (y-y^3) for this example.

 f[i] = 0.0;
 for( size_t j = 0; j < num_points; j++ ){
 f[i] += j*(j-1)*pow(y[j+num_points],j-2)*y[j]}
 f[i] *= std_noise;

if (i > 0){ for( size_t j = 0; j < num_points; j++ ){ f[i] += i*(y[j+num_points]-pow(y[j+num_points],3.0))*pow(y[j+num_points],i-1)*y[j];} }

 growth[i] = (y[i+num_points]-pow(y[i+num_points],3.0));

Return control to the calling routines.

 return;    
 }