CFD Online Discussion Forums - Help Parallelizing UDF

Hello members, this is my first time parallelizing an UDF and I'd like to know if there's somehting wrong as far as parallel syntax implementation. This code works properly in serial without all of the added parallel syntax. I should also mention the UDF compiles properly, but the solution diverges within the first time-step. Here is a brief description of the simulation setup and what I'm trying to accomplish:

Solver: 2D, laminar, incompressible, unsteady.
Computational Domain: Y-pipe with one inlet and two outlets.
Purpose of Model: To control the mass flow rate at the outlets using pressure defined by a UDF.
Inlet: unsteady mass-flux inlet with a somewhat sinusoidal profile resembling a heart beat.
Outlets: the UDF attempts to prescribe pressure at the outlets. This pressure is function of the inlet pressure and two other parameters, one of which is user defined with the purpose of controlling the fraction of total flow going out of the outlet.
Initialization: The simulation is initialized with a steady-state solution with equivalent inlet mass-flow rate.

Here is a brief description of the UDF:

-The function VelInlet specifies the unsteady mass flow rate.
-The define_adjust function gathers information such as pressure and mass flow rate from the boundaries. It also calculates the new values to be imposed.
-The define_profile functions impose the boundary conditions at the inlet(VelInlet mass flow rate) and outlets (average pressure calculated in define_adjust function). The pressure outlets are denominated RSA and LSA.

Here is the UDF

Code:

/* Udf for Inlet Velocity and BC prescription*/

 

#include "udf.h"

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

 

static real Q = 800.0*0.951*(1.0/(1000.0*1000.0*60.0))*1060.0;

static real pi = 3.1415926535898;

static real BPM = 130.0;

static real inlet_area = 0.0003125;

static real R_RSA = 10000.0, R_LSA = 40000.0;

 

int id_RSA,id_RCA,id_LSA,id_LCA,id_DA,id_RPA,id_LPA,id_RcorA,id_LcorA;

 

static real mflow, PressInlet, PressRSA, PressLSA, R_cfd_RSA, R_cfd_LSA, R_eq ;

static real mdot_RSA, mdot_LSA, time;

static real mflow_previous; /*For time step size criteria*/

static real sum1=0.0, sum2=0.0, sum3=0.0, AreaTot=0.0, AreaTot2=0.0, AreaTot3=0.0;

 

 

static real VelInlet(real t){

static real massflow;

static real sum = 0.0;

static real omega;

omega = pi*BPM;

static real A[8];

A[0] = 2.0*Q;

A[1] = -0.0235214;

A[2] = 0.0156809333333333;

A[3] = -0.0078404666666667;

A[4] = 0.0028510787878788;

A[5] = -0.0007127696969697;

A[6] = 0.0001096568764569;

A[7] = -0.0000078326340326;

static real B[8];

B[0] = 0.0;

B[1] = 0.0103416322216163;

B[2] = -0.013788842962155;

B[3] = 0.0103416322216163;

B[4] = -0.0050141247135109;

B[5] = 0.0015669139729722;

B[6] = -0.0002892764257795;

B[7] = 0.000024106368815;

 

int i;

sum = 0.0;

for (i = 1; i<8; i++){

sum += A[i]*cos(2.0*i*omega*(t/60.0))+B[i]*sin(2.0*i*omega*(t/60.0));

}

massflow = (A[0]/2.0 + sum)/inlet_area;

/* if (massflow=0.0){*/

/* To offset the flow from zero so that solver dosen't blow up,

offset to mflowrate of 0.001 do 0.001/0.0003125 = 3.2 */

/*massflow = 3.2;

}*/

 

return massflow;

}

 

 

DEFINE_DELTAT(mydeltat,d)

{

static real time_step;

if(fabs(mflow)<0.01){

time_step = 0.0002;

}

else if((fabs(mflow_previous) + fabs(mflow))< 1e-5){

time_step = 0.05;

}

else{

time_step = 0.001;

}

return time_step;

}

　

 

DEFINE_ADJUST(boundary_data,d)

{ 

time = CURRENT_TIME + 0.21;

mflow = 0.0;

mflow = VelInlet(time)*inlet_area;

mflow_previous = VelInlet(time-0.02)*inlet_area; /*For time step size criteria*/

 

sum1=0.0, sum2=0.0, sum3=0.0, AreaTot=0.0, AreaTot2=0.0, AreaTot3=0.0;

 

#if !RP_HOST

static real NV_VEC(A);

int zone_ID = 6;

face_t f;

Domain *domain;

domain = Get_Domain(1);

Thread *t = Lookup_Thread(domain, zone_ID);

 

begin_f_loop(f,t) 

if (PRINCIPAL_FACE_P(f,t)) 

{ 

F_AREA(A,f,t);

AreaTot += NV_MAG(A);

sum1 += NV_MAG(A)*F_P(f,t);

}

end_f_loop(f,t)

 

#if RP_NODE

sum1 = PRF_GRSUM1(sum1);

AreaTot = PRF_GRSUM1(AreaTot);

#endif /*RP_NODE*/

 

mdot_RSA = 0.0;

 

static real NV_VEC(B);

int ID_outRSA = 5;

id_RSA = ID_outRSA;

Thread *tRSA = Lookup_Thread(domain,ID_outRSA);

 

begin_f_loop(f,tRSA) 

if (PRINCIPAL_FACE_P(f, tRSA))

{

mdot_RSA = F_FLUX(f,tRSA)+ mdot_RSA;

}

end_f_loop(f, tRSA)

 

begin_f_loop(f,tRSA) 

if (PRINCIPAL_FACE_P(f, tRSA))

{ 

F_AREA(B,f,tRSA);

AreaTot2 += NV_MAG(B);

sum2 += NV_MAG(B)*F_P(f,tRSA);

}

end_f_loop(f,tRSA)

 

#if RP_NODE

mdot_RSA = PRF_GRSUM1(mdot_RSA);

sum2 = PRF_GRSUM1(sum2);

AreaTot2 = PRF_GRSUM1(AreaTot2);

#endif /*RP_NODE*/

 

mdot_LSA = 0.0;

 

static real NV_VEC(C);

int ID_outLSA = 4;

id_LSA = ID_outLSA;

Thread *tLSA = Lookup_Thread(domain,ID_outLSA);

 

begin_f_loop(f,tLSA) 

if (PRINCIPAL_FACE_P(f, tLSA))

{

mdot_LSA = F_FLUX(f,tLSA)+ mdot_LSA;

}

end_f_loop(f, tLSA)

 

begin_f_loop(f,tLSA) 

if (PRINCIPAL_FACE_P(f, tLSA)) 

{ 

F_AREA(C,f,tLSA);

AreaTot3 += NV_MAG(C);

sum3 += NV_MAG(C)*F_P(f,tLSA);

}

end_f_loop(f,tLSA)

 

#if RP_NODE

mdot_LSA = PRF_GRSUM1(mdot_LSA);

sum3 = PRF_GRSUM1(sum3);

AreaTot3 = PRF_GRSUM1(AreaTot3);

#endif /*RP_NODE*/

#endif /*!RP_HOST*/

　

node_to_host_real_6(sum1, sum2, sum3, AreaTot, AreaTot2, AreaTot3);

node_to_host_real_2(mdot_RSA, mdot_LSA);

 

#if !RP_NODE

PressInlet = sum1/AreaTot;

PressRSA = sum2/AreaTot2;

R_cfd_RSA = (PressInlet - PressRSA)/mdot_RSA;

PressLSA = sum3/AreaTot3;

R_cfd_LSA = (PressInlet - PressLSA)/mdot_LSA;

R_eq = 1.0/(1.0/(R_cfd_RSA + R_RSA) + 1.0/(R_cfd_LSA + R_LSA));

#endif /*RP_NODE*/

 

host_to_node_real_6(PressInlet, PressRSA, PressLSA, R_cfd_RSA, R_cfd_LSA, R_eq);

 

#if !RP_NODE

Message("Inlet Pressure: %g\n", PressInlet);

Message("RSA Pressure: %g\n", PressRSA);

Message("LSA Pressure: %g\n", PressLSA);

Message("Inlet Massflow: %g\n", mflow);

Message("Rcfd RSA: %g\n", R_cfd_RSA);

Message("Rcfd LSA: %g\n", R_cfd_LSA);

Message("R_eq: %g\n", R_eq);

#endif /*RP_NODE*/

}

　

 

DEFINE_PROFILE(massflux_inlet,thread,i)

{

 

#if !RP_HOST

int zone_ID = 6;

face_t f;

Domain *domain;

domain = Get_Domain(1);

Thread *t = Lookup_Thread(domain, zone_ID);

 

begin_f_loop(f,t) 

if (PRINCIPAL_FACE_P(f, t)) 

{

F_PROFILE(f,t,i) = mflow/inlet_area;

}

end_f_loop(f,t)

 

#endif /*!RP_HOST*/

}

　

DEFINE_PROFILE(RSA_Pressure, thread, i)

{

#if !RP_HOST

int ID_outRSA = 5; 

face_t f;

Domain *domain;

domain = Get_Domain(1);

 

id_RSA = ID_outRSA;

Thread *tRSA = Lookup_Thread(domain,ID_outRSA);

 

begin_f_loop(f,tRSA) 

if (PRINCIPAL_FACE_P(f, tRSA))

{

F_PROFILE(f, tRSA, i) = PressInlet - R_eq*mflow + R_RSA*mdot_RSA;

}

end_f_loop(f, tRSA)

 

#endif /*!RP_HOST*/

}

DEFINE_PROFILE(LSA_Pressure, thread, i)

{

 

#if !RP_HOST

int ID_outLSA = 4;

face_t f;

Domain *domain;

domain = Get_Domain(1);

 

id_LSA = ID_outLSA;

Thread *tLSA = Lookup_Thread(domain,ID_outLSA);

 

begin_f_loop(f,tLSA) 

if (PRINCIPAL_FACE_P(f, tLSA))

{

F_PROFILE(f, tLSA, i) = PressInlet - R_eq*mflow + R_LSA*mdot_LSA;

}

end_f_loop(f, tLSA) 

 

#endif /*!RP_HOST*/

}

-The function VelInlet specifies the unsteady mass flow rate.
-The define_adjust function gathers information such as pressure and mass flow rate from the boundaries. It also calculates the new values to be imposed.
-The define_profile functions impose the boundary conditions at the inlet(VelInlet mass flow rate) and outlets (average pressure calculated in define_adjust function).

This is the simplest case I was able create as proof of concept for a larger more involved model, so I apoligize if its a bit tedious.

Thanks very much for your time and attention.

AndresC.
Graduate Student,
University of Central Florida
MMAE