/***************************************************************************/
/* Implementation of                                                       */
/*                                                                         */
/*  (1) the unstready term for a user-defined scalar (1st-order)           */ 
/*                                                                         */
/*  (2)mixed boundary condition for a UDS (or multiple UDS's):             */
/*       q = hC ( phi - PHI_inf )                                          */
/***************************************************************************/

#include "udf.h"
#include "sg.h"  /* needed for the boundary and secondary gradient  macros */


/*=======================================================================*/
#define CP 1650.0  /* heat capacity for Maple in (J/kg K) */
#define HTC 20.0   /* heat transfer coefficient for the problem */
#define TINF 550   /* ambient temperature of the problem */
/*=======================================================================*/

/* Names of the user-defined scalar to be used */
enum
  {
    phi1,
    N_REQUIRED_UDS
  };

DEFINE_UDS_UNSTEADY(unst1stOrder, c, t, i, apu, su)
{
  /* if the unsteady term is different from the default term: d(rho*phi)
                                                              ----------
                                                                  dt
     this macro is used to specify the appropriate unstready term         */

  real volume, cp=CP, deltaTime=CURRENT_TIMESTEP;

  volume=C_VOLUME(c, t);

  /* the transient term is moved to the RHS of the equation and is split 
    into two parts---check the FVM algorithm for detail. First-order 
    backward differencing is implemented below:                           */

  *apu = -cp*volume/deltaTime;  
  *su  =  cp*volume*C_UDSI_M1(c, t, i)/deltaTime;

}


DEFINE_UDS_UNSTEADY(unst2ndOrder, c, t, i, apu, su)
{
  real volume, cp=CP, deltaTime=CURRENT_TIMESTEP;
  real temp;
  volume=C_VOLUME(c,t);

  /* printf("Time Step= %d\n", N_TIME); */

  if (N_TIME <  2)
    unst1stOrder(c, t, i, apu, su);
  else
    {
      /* 2nd-order backward differencing: 
	 c d phi/dt =  (1.5phi - 2phi^(n-1) + 0.5 phi^(n-2))*c*volume/detaT */
 
      temp=cp*volume;
      *apu = -1.5*temp/deltaTime;

      *su  = (2.0*C_UDSI_M1(c, t, i) - 0.5*C_UDSI_M2(c, t, i))*temp/deltaTime;

    }
}
DEFINE_PROFILE(scalarMixedBC, thread, nv)
{
  /* constants must be specified correctly for the mixed BC in the preamble section */
  real hC, PHI1_inf;

  /* ==================== */

  face_t f;
  
  real A[ND_ND], dG[ND_ND], dr0[ND_ND], es[ND_ND], dr, A_by_es;
  real Af;
  real beta0, gamma;
  real temp1, temp2;
  

  Thread *t0=thread->t0; 
  hC=HTC;
  PHI1_inf=TINF;


  begin_f_loop(f, thread)
    {

     /* identify the cell thread adjacent to the face thread f */
      cell_t c0 = F_C0(f, thread); 
 
      BOUNDARY_FACE_GEOMETRY(f, thread, A, dr, es, A_by_es, dr0);
      Af=NV_MAG(A);
      gamma=C_UDSI_DIFF(c0, t0, phi1);
 
      if (NULLP(T_STORAGE_R_NV(t0, SV_UDSI_G(phi1))))
	beta0=0;  /* if gradient is not allocated and stored yet, bypass
                     the following macro (it happens when saved
                     case/data files are being read during read-case-data */
      else
        BOUNDARY_SECONDARY_GRADIENT_SOURCE(beta0, SV_UDSI_G(phi1),dG, es, A_by_es, gamma);

  /* temporary variables used in the profile expression */
      temp1=gamma*A_by_es/dr;
      temp2=hC*Af;
      
      
      F_PROFILE(f, thread, nv)
	= (temp1*C_UDSI(c0, t0, phi1)- beta0 + temp2*PHI1_inf)/(temp2 + temp1);
    }
  end_f_loop(f, thread)
  
}
