Question regarding flowrates

I have implemented open boundary condition according to Zou and He at
the inlet and at the outlet of my system;
I have imposed a poiseuille velocity profile at the inlet, v.
How can I compute mass flow rates, Qin and Qout?

the mass flow rate is the integral (sum) of density times velocity. You could compute it like this:

double mass_flow_x = 0; // mass flow in x-direction
X = 0 // inlet or
// X = X_max // outlet

for(Y = 0; Y < Y_max; ++Y) {
for(Z = 0; Z < Z_max; ++Z) {
mass_flow += density(X, Y, Z) * velocity_x(X, Y, Z)

By the way, if mass_flow_x at inlet and outlet are not the same, the total mass in the system will change.