[Neb]

Hi everyone,
i have to write a new BC in openfoam and i am using codedFixedValue. I have already asked for help in other posts and tried a lot, but the situation does not improve.
I've been stuck for a long time and I don't know how to proceed.
Can anyone tell me how to read the outlet flow rate on the outlet at a certain moment?


const scalar &t = this ->db().time().value(); //time
const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi"); //flow rate
const fvsPatchField<scalar>& phip=patch().patchField<surfaceScalarField,scalar> (phi);
float i= t;
const scalar dt=0.001;

for (i=0; i<=1; i+=dt);
{ const scalar C=1.16e-05;
const scalar R1=1.64e+03;
const scalar R2=1.79e+05;
const scalar p_0=7.55;
const scalar Q_0=1.03e-4;
static scalar integrand (0);
integrand +=(exp(i/(R2*C))*gSum(phip)*dt);
operator ==(R1*gSum(phip)+exp(-i/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));
}
endl;


In this code I have the following problems:
1) I run the simulation for two seconds, but the time must remain in the first second. Is the for loop well written?

2) The integrand function increases the outlet flow rate for both simulation seconds, but I want it finishes at the first second, add the flow rate from t = 0 and not to t = 1 + dt.

Please, someone help me, I don't know how to solve anymore.

[RGS]

Originally Posted by Neb

Hi everyone,
i have to write a new BC in openfoam and i am using codedFixedValue. I have already asked for help in other posts and tried a lot, but the situation does not improve.
I've been stuck for a long time and I don't know how to proceed.
Can anyone tell me how to read the outlet flow rate on the outlet at a certain moment?

I cannot unfortunately help you with the code that you posted because I am still getting familiar with how to code a custom boundary condition. But shouldn't reading the flow rate at the outlet be similar to reading the value at any other point in the mesh? Perhaps you might have to interpolate from he cell center to the surface to get the exact value at the patch. I think I can find the code for that with a bit of searching.

[Neb]

Hi RGS,

Thank you for your answer. I need to read the flow rate in a given instant. With gSum (phip) I know the range on the patch outlet, I just have to be able to recall the one calculated in a certain t. But since I'm stuck on this point, I'm writing the code in another way, I hope it works.

[crubio.abujas]

Hello again Neb!

I remember we've talked about this integral approach, but I have the sense that you're a little bit confused about what you're code is really doing. Let me try to explain it bit by bit so you can better understand it.

First, you have to know when this specific code is executed. The solver is instructed to first set the coefficients for the interior cell equation and afterwards update the coefficients using the boundaries definitions, so this code will be executed each time the solvers instruct it to do it. Depending on the solver you're using this may happen once per iteration or several times along the solving procedure.

I am not totally sure what kind of application you're using this code for, but I'm going to guess is not a steady-state problem with a pseudo-time or any other fancy approach. In that case, the time management is on the solver-side. The job of the boundary condition is, when asked by the solver, to provide a valid set of patch values to set the equations. This values can depend on time, but this time is provided by the solver, not the other way around.

I've added a commented version of the code in case it help you better to understand it:

const scalar &t = this ->db().time().value();

const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi"); 
const fvsPatchField<scalar>& phip=patch().patchField<surfaceScalarField,scalar> (phi);

// This line is assigning to i the value of the actual t. 
// It it also doing an implicit transformation as, under the hood, scalar =
// double, so you're loosing precision
float i= t; 
const scalar dt=0.001;

// On this loop you're rewriting the value of i, now i and t are totally
// independent. Changing the value of i will not change the value of t in any
// mean. So This can be interpreted as a "pseudo-time" from 0 to 1.
for (i=0; i<=1; i+=dt); { 

    // Why are you defining all these constants inside the loop?
    // The scope of these is the entire function, so its a better approach to
    // set them once outside the loop.
    const scalar C=1.16e-05;
    const scalar R1=1.64e+03;
    const scalar R2=1.79e+05;
    const scalar p_0=7.55;
    const scalar Q_0=1.03e-4;

    // The idea of setting a scalar called integrand is to tell OF to keep the
    // variable value even when it goes outside the BC function. But in this
    // case I think you are not using this feature to pass information across
    // diferent instants as i is being reset on each call
    static scalar integrand (0);
    integrand +=(exp(i/(R2*C))*gSum(phip)*dt);

    // Again, it have no sense to set the values of current patch inside the
    // loop, as they are not going to be used until you return the control to
    // the solver. So you're just overwriting the value on each new loop and
    // the solver uses the last instant iteration.
    operator ==(R1*gSum(phip)+exp(-i/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));
}

Concerning your questions:

In this code I have the following problems:
1) I run the simulation for two seconds, but the time must remain in the first second. Is the for loop well written?

As I've explained you. You cannot alter the time from inside a boundary condition. That has to be done from the solver controls.

2) The integrand function increases the outlet flow rate for both simulation seconds, but I want it finishes at the first second, add the flow rate from t = 0 and not to t = 1 + dt.

I don't think I've comprehend this question. I think it may be related with the willing of finish the simulation when an specific time is reached, but I dont get what do you mean with "both simulation seconds".

I hope that explanation helps you reframe your approach. What do you want to achieve with such code?

[Neb]

hello again ,
thank you so much for the explanation. You are very kind. It is difficult to explain what I want to do. The first code doesn't return the pressure values ​​that I expect, so I had to make some changes that work, but they're still not perfect. I try to explain myself again. Work with pimple, unsteady flow. In the code I added the for loop because I want the time to remain between 0 and 1 (and it worked), while the simulation time will run between 0 and 2. Same thing for the sum: it must add values ​​between 0 and 1, but it must start again between 1 and 2. It must not have memory of the flow between 0 and 1 second. I thought of a way to solve this too, but I have yet to check if it works.
I will read carefully what you have explained to me. Thanks again for your time.

[crubio.abujas]

The problem with the loop approach is that you're ALWAYS evaluating the integrand between 0-1s. Even when you are at t=0.1 or t=0.2 and there is still no information on the reaming points. Additionally, when you evaluate the boundary condition you only have access to the current phi values (or previous ones if you've stored them, but you have to recall them explicitly). So the evaluation of the integral has been done using wrong values.

It seems that you are trying to apply some kind of outlet saturation. Maybe you're trying to model some kind of vessel and how its pressure evolve along the filling up process. You may be interested in limiting the integral part, as when a given point is achieved the integrand may become huge and cause problem, but the overall effect is negligible because of the exponential factor multiplying it. If that is the case, you may simply add an If statement to stop calculating when an specific time is reached.

if (t < T_REF) {
    integrand +=(exp(t/(R2*C))*gSum(phip)*dt);
    operator ==(R1*gSum(phip)+exp(-i/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));
}

With this code you evaluate the equation while t is lower than a reference value, when t is higher it does no further calculation and keeps the last calculated value.

But maybe this is not the expected behavior. You mentioned that you want the conditions to be restarted at t=1s. If that is the case, you can just check if the condition is reached and restart properly:

if (t > T_REF)  {
     integrand = 0;
     // other resetting options
}
integrand +=(exp(t/(R2*C))*gSum(phip)*dt);
operator ==(R1*gSum(phip)+exp(-i/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));

This will create a discontinuous jump on the pressure, but maybe this is what you expects.


Let me know if that is close to your intentions!

[Neb]

I carefully read your explanations in the first post and now some things are clearer. Yes, you well hypothesized, I am modeling a vassel and I need this pressure correction to model the missing downstream arteries that offer resistance to the flow in my artery.

I knew that the time is dictated by the solver, but I need that in the BC function it always remains between 0 and 1 second (real heart cycle), regardless of the endTime set in the solver (which are 2 seconds = two cardiac cycles).

Exactly, the integral defined initially creates enormous pressure values. So I need to create a discontinuity at 1 second, without adding the calculated flow rate from 0 to 0.9 second.

The correction I had hypothesized for the integral is exactly the IF correction instruction that you suggested to me.
I thought such a thing: if t = 1 + dt computes the integral for this value and starts adding again. But I can't translate it well.

[crubio.abujas]

What an interesting application! If you're planning to run other duration simulation you may consider to use modulus operator on a pseudo-time:

scalar tau = t % T_REF;

// Reset condition
if (tau < dt) {
    integrand = 0;
}
integrand +=(exp(tau/(R2*C))*gSum(phip)*dt);
operator ==(R1*gSum(phip)+exp(-tau/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));

The modulus operator (%) will return a value between 0 and T_REF (1 in your case), so for example if you get a t value of 2.3s it will return 0.3s. The integrand is the only parameter that store information of previous steps so when a new cycle starts (tau < dt) it is set to 0 again. It is better to use ranges instead of precise values as "tau == 0" because you may have numerical errors or different timesteps that give you some unexpected behaviour. Is more general this way. With this approach you may have endTimes of 3, 4 or whatever number of secods you may like and the boundary condition will only operate among 0-1 tau values.

I was wondering about p_0. Is this a constant value or is referred to the initial step? I mean, shall you update p_0 on the next cycle or you just set it on the beginning and change it no more?

[Neb]

yes, It is a very particular application.

Great!! Exactly, a generic approach is perfect. This way I can do all the heart cycles I want.

The value of p_0 is delicate. By doing the sensitivity study of the parameters of the equation, I realized that it greatly changes the range of pressure values. So I have to evaluate it carefully. Currently I enter an average value and do not change it, but the behavior must be well evaluated by running the simulation.

[crubio.abujas]

Good to know it suits your needs! Let me know any further advance you have with the model!

[Neb]

Sorry for the delay in response, I will let you know if I get the desired pressure values! Thanks again

[Neb]

By entering the new code It returns the following error: "invalid operands of type" scalar const {aka const double "and" int "in binary" operator% "
it refers to the row:

scalar tau = t % T_REF

[crubio.abujas]

Sorry Neb, that was my fault. Although "%" is the modulus operator, it is only valid for integers variables (not double or scalar). I've checked the concept into python and assumed that it will be the same in C++.

Luckily it has an easy solution. Just replace it by:

scalar tau = std::fmod(t, T_REF);

It worked for me (this time I've tried in OF ), hopefully it will work for you!

[Neb]

yes it works now!
No problem, you helped me a lot.
I tried to solve it without disturbing you again, looking for an alternative method, but obviously I'm too inexperienced with this language to do it.: o

[crubio.abujas]

Don't worry. This is helping me to better understand OF as well. Nevertheless is always a good things that you find a way to solve problems, so I encourage you to persist in that path. Of course, if you get stuck at some point there is no shame in asking for help!

Good luck with your model

[Neb]

Hello,

unfortunately there is still something wrong

I realized that the integral causes enormous pressure values yet. Eliminating gSum (phip), the final pressure values ​​are better.
I hypothesized that the problem lies in the + = increment. That operation takes all the flow values ​​at each iteration of the specific t and adds them together. Instead it has to add up the final flow of each time t, t+dt etc.. not the flow rate at each iteration. I don't know if I managed to make myself understood and if it is true that this happens...

This is the part that still doesn't work well:

integrand +=(exp(tau/(R2*C))*gSum(phip)*dt);

To correct how can I introduce the flow rate at the time I want instead of gSum(phip)? Is there any way?

[crubio.abujas]

Well yes, that can definitely be the case. You can add a variable following the time and only update the integrand then the current_time and t are different. That way you can leave the increment only for real timeseps. I haven't tested this code, so maybe it need some debugging, but you can try with this approach. Let's hope it works!

const scalar &t = this ->db().time().value();

const surfaceScalarField& phi = db().lookupObject<surfaceScalarField>("phi"); 
const fvsPatchField<scalar>& phip =
    patch().patchField<surfaceScalarField,scalar> (phi);
const fvsPatchField<scalar>& phip_last =
    patch().patchField<surfaceScalarField,scalar> (phi.oldTime());

// Constants
const scalar C=1.16e-05;
const scalar R1=1.64e+03;
const scalar R2=1.79e+05;
const scalar p_0=7.55;
const scalar Q_0=1.03e-4;

// Static variables
static scalar integrand (0);
static scalar current_time (0);

scalar tau = std::fmod(t, T_REF);

 // Check if the timestep in new
if (current_time < t) {

    // Update the stored timestemp
    current_time = t;

    // Update the integrand with the previous value
    integrand += (exp(tau/(R2*C))*gSum(phip_last)*dt);

    // Check for reset condition
    if (tau < dt) {
        integrand = 0;
    }
}

operator ==(R1*gSum(phip)+exp(-tau/(R2*C))*((p_0-R1*Q_0)+(1/C)*integrand));

[Neb]

Hello, first of all thanks again for the suggestion.

I tried the code but the loop is not used. the integral works with gSum (phip_last), the first cardiac cycle from 0 to 1 works well, but the second cycle with the time from 1 to 2, no. I don't understand what the code is doing now since it no longer reads the time corrections that were made previously.

[crubio.abujas]

Originally Posted by Neb

Hello, first of all thanks again for the suggestion.

I tried the code but the loop is not used. the integral works with gSum (phip_last), the first cardiac cycle from 0 to 1 works well, but the second cycle with the time from 1 to 2, no. I don't understand what the code is doing now since it no longer reads the time corrections that were made previously.

What kind of error are you finding?
Maybe the reset condition is not working properly. Try it wirh:

If (tau <= dt) {
     integrand = 0 ;

}  

Info << "tau = " << tau <<  " - dt = " << dt << endl;

The second part is to print the value of tau to check the condition.

[Neb]

The code compiles, but says that the loop is unused.
In fact, removing the loop, only this command line is read:
integrand + = exp(t/R2*C)*gSum(phip_last)*dt.

Any corrections that I am going to make according to time do not consider it.

Now I try your new suggestion.