[Morice]

Hey guys,

The UDF is an energy and water vapour source UDF and I keep getting a floating-point error and divergence on temperature and water vapour when I try to run the simulation.

I have done the math in C++ and Matlab and the math is correct.

Maybe there is an issue with the thermodynamics that am not seeing.

Honestly, I have no idea what the issue could be.

Thanks so much, much appreciated.

 #include "udf.h"
#include "math.h"

#define ra 271           /*aerodynamic resistance*/
#define lad .3033        /*leaf area demsity*/
#define kc  .95        /*extinction coefficient of the canopy*/
#define hgt 1.5         /*height of the canopy*/
#define Cp  1006.43     /*specific heat of dry air*/
#define rho 1.225       /*density of dry air*/
#define lai .24271      /*leaf area index*/
#define c3  18.6
#define c4  197.5
#define c5  .31
#define c6  1.2e-6
#define lamda 2.450e6       /*latent heat of vaporisation*/
#define gamma 67.381125     /*psychrometric constant*/
#define ks    0.7           /*coefficient of shading screen*/



/*VARIABLES THAT NEED TO BE MODIFIED BEFORE EVERY SIMULATION*/

#define rgo  43.28        /*global radiation above the canopy*/

DEFINE_INIT(initial_values, domain)
{
	cell_t c;
	Thread* t;
	thread_loop_c(t, domain)
	{
		begin_c_loop(c, t)
		{
			C_UDMI(c, t, 0) = 288.15;           /*Ta*/
			C_UDMI(c, t, 1) = 1699.81;         /*Psat*/
			C_UDMI(c, t, 2) = 289.15;          /*Temperature of leaf in cropzone1*/
			C_UDMI(c, t, 3) = 60.497;          /*rs1*/
			C_UDMI(c, t, 4) = 331.495;         /*rt1*/
			C_UDMI(c, t, 5) = 6.19434;         /*ET1, Qlat1*/
			C_UDMI(c, t, 6) = 292.226;         /*Tlr1; Real leaf temp cropzone1*/
			C_UDMI(c, t, 7) = 7.66834e-7;      /*cropzone1 water vapour mass fractions source*/
			C_UDMI(c, t, 8) = 289.15;          /*Temperature of leaf in cropzone2*/
			C_UDMI(c, t, 9) = 60.4947;         /*rs2*/
			C_UDMI(c, t, 10) = 331.495;         /*rt2*/
			C_UDMI(c, t, 11) = 7.66834e-7;      /*cropzone2 water vapour mass fractions source*/
			C_UDMI(c, t, 12) = 292.226;         /*Tlrl; Real leaf temp crop zone2*/
			C_UDMI(c, t, 13) = 6.19434;         /*ET2, Qlat2*/
			
		}
		end_c_loop(c, t)
	}
}

DEFINE_ADJUST(ta_calc, d)
{
	Thread* t;
	cell_t c;
	thread_loop_c(t, d)
	{
		if(THREAD_ID(t)==7)
		begin_c_loop(c,t)
		{
			double Ta,Psat;
			
			Ta += C_T(c, t) * C_VOLUME(c,t);
			
			C_UDMI(c, t, 0) = Ta;
			
			Psat = 610.78 * exp((17.25 * (Ta - 273.15)) / (237.8 + (Ta - 273.15)));
			
			C_UDMI(c, t, 1) = Psat;
		}
		end_c_loop(c,t)
	}
}


DEFINE_SOURCE(cropzone1_h20_source, c, t, dS, eqn)
{
	double a, rs1, rt1, hcv1, Tl1, Tl1r, Plsat1, vpd1, ET1, Sw1, St1;

	real con, source;

	real x[ND_ND];

	C_CENTROID(x, c, t);

	con = rgo * ks * exp(-kc * lad * (((hgt - x[1]) / hgt)));

	Tl1 += C_T(c, t) * C_VOLUME(c, t);

	C_UDMI(c, t, 2) = Tl1;

	Plsat1 = 610.78 * exp((17.25 * (Tl1 - 273.15)) / (237.8 + (Tl1 - 273.15)));

	vpd1 = Plsat1 - C_UDMI(c, t, 1);

	a = con / (2 * lai);

	rs1 = c3 * ((a + c4) / (a + c5)) * (1 + (c6 * vpd1 * vpd1));

	C_UDMI(c, t, 3) = rs1;

	rt1 = ra + rs1;

	C_UDMI(c, t, 4) = rt1;

	ET1 = ((rho * Cp) / gamma) * (vpd1 / rt1);

	C_UDMI(c, t, 5) = ET1;

	hcv1 = (rho * Cp) / ra;

	Tl1r = (con - ET1) / ((2 * hcv1) + C_UDMI(c, t, 0));

	C_UDMI(c, t, 6) = Tl1r;

	source = (rho * Cp * lad * vpd1) / (gamma * lamda * rt1);

	C_UDMI(c, t, 7) = source;

	dS[eqn] = (rho * Cp * lad) / (gamma * lamda * rt1);

	return source;
}


DEFINE_SOURCE(cropzone1_energy_source, c, t, dS, eqn)
{
	real source;

	source = - C_UDMI(c, t, 7) * lad;

	dS[eqn] = -(rho * Cp * lad) / (gamma * C_UDMI(c, t, 4));

	return source;
}


DEFINE_SOURCE(cropzone2_h20_source, c, t, dS, eqn)
{
	double b, rs2, rt2, hcv2, Tl2, Tl2r, Plsat2, vpd2, ET2, Sw2, St2;

	real x[ND_ND];

	real source, con;

	C_CENTROID(x, c, t);

	con = rgo * ks * exp(-kc * lad * (((hgt - x[1]) / hgt)));

	Tl2 += C_T(c, t) * C_VOLUME(c, t);

	C_UDMI(c, t, 8) = Tl2;

	Plsat2 = 610.78 * exp((17.25 * (Tl2 - 273.15)) / (237.8 + (Tl2 - 273.15)));

	vpd2 = Plsat2 - C_UDMI(c, t, 1);

	b = con / (2 * lai);

	rs2 = c3 * ((b + c4) / (b + c5)) * (1 + (c6 * vpd2 * vpd2));

	C_UDMI(c, t, 9) = rs2;

	rt2 = ra + rs2;

	C_UDMI(c, t, 10) = rt2;

	ET2 = ((rho * Cp) / gamma) * (vpd2 / rt2);

	C_UDMI(c, t, 11) = ET2;

	hcv2 = (rho * Cp) / ra;

	Tl2r = (con - ET2) / ((2 * hcv2) + C_UDMI(c, t, 0));

	C_UDMI(c, t, 12) = Tl2r;

	source = (rho * Cp * lad * vpd2) / (gamma * lamda * rt2);

	C_UDMI(c, t, 13) = source;

	dS[eqn] = (rho * Cp * lad) / (gamma * lamda * rt2);

	return source;
}

DEFINE_SOURCE(cropzone2_energy_source, c, t, dS, eqn)
{
	real source;

	source = - C_UDMI(c, t, 13) * lad;

	dS[eqn] = -(rho * Cp * lad) / (gamma * C_UDMI(c, t, 10));

	return source;
}