Re: [DuMux] Question on sourceAtPos function

2021-04-08 Thread Timo Koch


> On 8. Apr 2021, at 13:42, Mahmoud Atef Mahmoud Mohamed Aboelseoud S277151 
>  wrote:
> 
> Dear Timo,
> 
> Thanks a lot for your reply. I didn't use the pointSources because I didn't 
> know how to assign the temperature inside such function for the injection 
> well because I want to impose a rate and also a temperature.

Hi Mahmoud,

for the pointSource for the OnePNI model you need to specify mass (kg/s) and 
energy source (J/s) (Indices::conti0EqIdx (=0) and Indices::energyEqIdx (=1)).
“NumEqVector" has size 2 because you have two equations.
The energy flow can be computed as massFlow*enthalpy(p,T) (neglecting heat 
conduction which is usually unimportant here). So the temperature enters via 
the enthalpy.
You can evaluate the correct enthalpy for your fluid system with 

using FluidSystem = GetPropType;
const auto enthalpy = FluidSystem::liquidEnthalpy(T, p);

for a temperature T (in K) and pressure p (in Pa). 

In your examples below you build a temporary object of type FluidState and set 
the temperature on this object. However this fluid state object is used 
nowhere, so it does nothing.
Consequently you seem to inject mass without accounting for the energy the 
injected fluid carries.

Hope this helps
Timo

> I'm using the OnePNI model in my problem. I will try to be more specific this 
> time as I'm afraid I could be making some mistake that I'm not aware of. To 
> define both the injection and production wells inside the sourceAtPos 
> function, I'm writing the following:
> 
> 
> NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
>{
> 
> NumEqVector values;
> using FluidState = GetPropType;
>FluidState fs;
> 
>// we define one source term and one sink term
>//Injection well (Source Term)
>if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) && 
> (globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && (globalPos[0] > 35 
> - eps_ && globalPos[0]  < 35 + eps_))
>{
>fs.setTemperature(293);
>values[contiWEqIdx] = 1e-1; // kg/(s*m^3)
>}
> 
> //Production well (Sink Term)
>else if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) && 
> (globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && (globalPos[0] > 65 
> - eps_ && globalPos[0]  < 65 + eps_))
>{
>values[contiWEqIdx] = -1e-1;  // kg/(s*m^3)
>}
> else {
> values[contiWEqIdx]=0;
> 
> }
> 
> return values;
> }
> 
> where contiWEqIdx is defined in an enum as contiWEqIdx = Indices::conti0EqIdx
> 
> and my model dimensions and cells are as follows:
> [Grid]
> Positions0 = 0 10 30 70 90 100
> Positions1 = 0 10 30 70 90 100
> Positions2 = 0 10 40 50
> Cells0 = 1 2 4 2 1
> Cells1 = 1 2 4 2 1
> Cells2 = 1 3 1
> 
> where there are upper and lower boundary layers with very low permeability 
> (order of e-25 m^2) each 10 m thick. I'm not using any refinement for the 
> moment around the wells. Is there anything wrong I'm making here in the above 
> function and which causes me to have a very limited range of mass rate below 
> which the simulation proceeds very quickly and above which I have really high 
> pressures and the simulator aborts ??!
> 
> In fact I want to impose a rate of 100 m^3/day in each well as it is a common 
> liquid rate and that's why I do the calculations mentioned in the previous 
> e-mail. but in order to implement what you proposed Timo in your response 
> which is using "elemVolVars[scv].density(phaseIdx)" for the density and 
> "scv.volume()" for the volume, shouldn't I divide the volumetric rate of the 
> well which is 100 m^3/day by (3 cells x 8 scv) before dividing by the volume 
> and multiplying by the density ? Anyway, I tried to follow what you proposed 
> by writing the following for the injection well for example
> 
> NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
>{
> SubControlVolume scv;
> ElementVolumeVariables elemVolVars;
> NumEqVector values;
> using FluidState = GetPropType;
>FluidState fs;
> 
>// we define one source term and one sink term
>//Injection well (Source Term)
>if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) && 
> (globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && (globalPos[0] > 35 
> - eps_ && globalPos[0]  < 35 + eps_))
>{
>fs.setTemperature(293);
>values[contiWEqIdx] = 
> (1.157e-3*elemVolVars[scv].density(LiquidIdx))/(24*scv.volume()); // 
> kg/(s*m^3)
>}
> 
> where LiquidIdx is defined in an enum as LiquidIdx = FluidSystem::comp0Idx
> 
> but I'm getting the following error:
> /home/mahmoud/dumux/dumux/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh:288:7:
>  note:   candidate expects 1 argument, 0 provided
> 
> 
> I would really appreciate it if you let me know how to correctly accurate for 
> the 100 m^3/day volumetric rate either using the 

Re: [DuMux] Question on sourceAtPos function

2021-04-08 Thread Mahmoud Atef Mahmoud Mohamed Aboelseoud S277151

Dear Timo,

Thanks a lot for your reply. I didn't use the pointSources because I 
didn't know how to assign the temperature inside such function for the 
injection well because I want to impose a rate and also a temperature. 
I'm using the OnePNI model in my problem. I will try to be more specific 
this time as I'm afraid I could be making some mistake that I'm not 
aware of. To define both the injection and production wells inside the 
sourceAtPos function, I'm writing the following:



NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
{

 NumEqVector values;
 using FluidState = GetPropTypeProperties::FluidState>;

FluidState fs;

// we define one source term and one sink term
//Injection well (Source Term)
if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) && 
(globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && (globalPos[0] 
> 35 - eps_ && globalPos[0]  < 35 + eps_))

{
fs.setTemperature(293);
values[contiWEqIdx] = 1e-1; // kg/(s*m^3)
}

 //Production well (Sink Term)
else if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) 
&& (globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && 
(globalPos[0] > 65 - eps_ && globalPos[0]  < 65 + eps_))

{
values[contiWEqIdx] = -1e-1;  // kg/(s*m^3)
}
 else {
 values[contiWEqIdx]=0;

 }

 return values;
 }

where contiWEqIdx is defined in an enum as contiWEqIdx = 
Indices::conti0EqIdx


and my model dimensions and cells are as follows:
[Grid]
Positions0 = 0 10 30 70 90 100
Positions1 = 0 10 30 70 90 100
Positions2 = 0 10 40 50
Cells0 = 1 2 4 2 1
Cells1 = 1 2 4 2 1
Cells2 = 1 3 1

where there are upper and lower boundary layers with very low 
permeability (order of e-25 m^2) each 10 m thick. I'm not using any 
refinement for the moment around the wells. Is there anything wrong I'm 
making here in the above function and which causes me to have a very 
limited range of mass rate below which the simulation proceeds very 
quickly and above which I have really high pressures and the simulator 
aborts ??!


In fact I want to impose a rate of 100 m^3/day in each well as it is a 
common liquid rate and that's why I do the calculations mentioned in the 
previous e-mail. but in order to implement what you proposed Timo in 
your response which is using "elemVolVars[scv].density(phaseIdx)" for 
the density and "scv.volume()" for the volume, shouldn't I divide the 
volumetric rate of the well which is 100 m^3/day by (3 cells x 8 scv) 
before dividing by the volume and multiplying by the density ? Anyway, I 
tried to follow what you proposed by writing the following for the 
injection well for example


NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
{
 SubControlVolume scv;
 ElementVolumeVariables elemVolVars;
 NumEqVector values;
 using FluidState = GetPropTypeProperties::FluidState>;

FluidState fs;

// we define one source term and one sink term
//Injection well (Source Term)
if ((globalPos[2] > 15 - eps_ && globalPos[2] < 35 + eps_) && 
(globalPos[1] > 45 - eps_ && globalPos[1]  < 45 + eps_) && (globalPos[0] 
> 35 - eps_ && globalPos[0]  < 35 + eps_))

{
fs.setTemperature(293);
values[contiWEqIdx] = 
(1.157e-3*elemVolVars[scv].density(LiquidIdx))/(24*scv.volume()); // 
kg/(s*m^3)

}

where LiquidIdx is defined in an enum as LiquidIdx = 
FluidSystem::comp0Idx


but I'm getting the following error:
/home/mahmoud/dumux/dumux/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh:288:7: 
note:   candidate expects 1 argument, 0 provided



I would really appreciate it if you let me know how to correctly 
accurate for the 100 m^3/day volumetric rate either using the 
sourceAtPos function with elaborating the mistakes or using the 
pointSource function showing how to impose the 293 kelvin temperature 
for the injection well.
I'm extremely sorry for the long e-mail. I just wanted to clarify my 
problem a bit more. Thanks in advance !



Best Regards,

Mahmoud Atef













Il 2021-04-07 20:27 Timo Koch ha scritto:

Dear Mahmoud,

what you are doing seems correct to me. You can get the relevant cell
volume with “scv.volume()” so you don’t have to compute that
yourself.
You can also get the density with
"elemVolVars[scv].density(phaseIdx)". But that shouldn’t change
anything if you already used the correct volume in your computation.

You can also use pointSources (see e.g
https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/test/porousmediumflow/1p/pointsources/timeindependent/problem.hh)
which is simpler for the case you describe.
There you just specify the position of the well(s) and the rate(s)
(this time directly in kg/s without the need