Cannot run the simulator

[nurlanj]

Hello,

sorry to bother you, I have read your paper and I noticed that in simulator you used different boundary conditions for TEO setup compared to ones that are in the paper. I assume you used grid-centered blocks with two ghost cells as boundaries. Also I assume you used harmonic means of mobilities with no Upwinding. I tried to use the boundaries you specified in the paper in my simulator but it didn’t converge even for the first time step. I would appreciate if you can explain me which boundary conditions to use.

In addition I tried to run your code but there is some kind of mistake which doesn’t allow me run it. Please could you check your uploaded files.

Thank you.

[jcabbasi]

Hello,

Thank you for reaching out. May I know what the problem is exactly (with my code I mean)?

Especially for the OEO case, I remember the boundaries are very critical and sensitive to implementation. To be honest, the code is a little messy; however, you can find the applied boundaries in the file 'FlowFncs.py'; function 'buildIt'. Can you read it properly?

BTW, I uploaded a new branch. Can you also check the new upload? 'Cent_Sim_Kumar2014 - 2.py' This file particularly.

[nurlanj]

Thanks for your reply and upload!

I will try to run your code again later this week.

In no way I am an expert in numerical simulations, so my questions may sound basic. If you look at your paper, eq. 36 for TEO boundaries and Python code (FlowFncs line 128 and 129) do not match a little bit, here I assume that 0 index in the paper is for innermost ghost cell. For me it looks like in your code you set this condition also to grid #1.

Another question is as you mention in eq. 35, Lambda_w (1/2)=0 (i.e. at the boundary), and theoretically it is only possible if both grids #0 and #1 are at Smin? The same is for outermost end of core, as non-wetting phase is immobile here, grid #N and #N+1 should be at Smax, as non-wetting phase mobility is zero at the boundary. Would it be possible to force set saturations in these cells to the abovementioned values and get right boundaries? Or I missed something?

[jcabbasi]

My Pleasure :) .

Regarding the 1st question: Cell [0] is the innermost cell (close to the center), and out of porous media. So, Pc=0 (infinite pore radius). Cell [1] can be assumed to be infinitely close to cell [0], but inside the porous media, so Pc is not zero. Line 128 addresses Eq. 36 when Pc=0. While Line 129 addresses Eq.36 when Pc is not zero. This is the difference between them.

Regarding question 2, I hope I understood it well. But at the innermost cell, where the water-saturated core is exposed to the nonwetting phase, no physical force allows water to be expelled from this side. However, we can not (to have better numerical performance) set the saturation cell [1] to be Smin. So, we should manually set the transmissibility to zero (Line 161). I am not sure about my answer, but I guess if you force-set the saturations, it may create some material balance problems. Setting the transmissibilities seems to be more physical than setting the saturations, from my perspective. However, you can quickly try both to see which one works better [I may not remember the details of my work correctly, since it is related to two years ago :) ).

[nurlanj]

Thank you very much for your reply) Now your logic is clear to me.

I installed all required modules and tried to run your code. There was an error in Visualization class, so I assumed that file name was different and renamed "dataVisualization_1" to "visualizationClass".

After this I got this error:

= RESTART: C:\Users\Nurlan\Downloads\centrifuge_public-Full_Code\centrifuge_public-Full_Code\Cent_Sim_Kumar2014.py Traceback (most recent call last): File "C:\Users\Nurlan\Downloads\centrifuge_public-Full_Code\centrifuge_public-Full_Code\Cent_Sim_Kumar2014.py", line 9, in import CentSimCore as simcore File "C:\Users\Nurlan\Downloads\centrifuge_public-Full_Code\centrifuge_public-Full_Code\CentSimCore.py", line 15, in from visualizationClass import * File "C:\Users\Nurlan\Downloads\centrifuge_public-Full_Code\centrifuge_public-Full_Code\visualizationClass.py", line 9, in optData = pd.read_pickle("./optData.pkl") File "C:\Users\Nurlan\AppData\Local\Programs\Python\Python311\Lib\site-packages\pandas\io\pickle.py", line 179, in read_pickle with get_handle( File "C:\Users\Nurlan\AppData\Local\Programs\Python\Python311\Lib\site-packages\pandas\io\common.py", line 868, in get_handle handle = open(handle, ioargs.mode) FileNotFoundError: [Errno 2] No such file or directory: './optData.pkl'

I assume there is a problem with optData.pkl missing in the directory.

How can I generate this file? Or there is some mistakes in the code? Sorry for so many questions) Hope you can help me.

[nurlanj]

Hello,

There is one more question regarding the boundaries:

For TEO case in the first grid you specified this boundary: PX[1]=self.PnwEqDuetoRotationMAX(omega)-rf.pc(SwX[1])

As the first grid cell is always updated during newton iterations, this boundary also changes, or you force set it to the specified values?

[jcabbasi]

Hello,

sorry for late reply. I added new files, please check the visualization class uploaded.

Regarding your question, since here the pressure is a function of Sw, it is updated in each iteration. However, since here the cell is in touch with the boundary, its Pressure value should be manually set at each iteration, as you see here.

[nurlanj]

Hello, thank you for uploaded new files. It was very helpful.

I used boundary conditions you mentioned earlier and now I can match theoretical saturation profile with numerical pretty close, but still there is some difference of about 1-2 saturation units on average. But when running your code it matches perfectly.

Could this be due to grid size that I use? My grid sizes are all the same for real and ghost cells. Also for calculations of terms in jacobian matrix, I mean partial derivatives, I use central difference, and for dx I use fixed value of about 1E-10. Could this also be a problem here?

Thank you in advance.

[jcabbasi]

Hi, I am happy it finally worked :).

Yes, one of the reasons can be the large grid sizes, especially at points close to the BCs (and small time steps at the beginning times). You can use large grid numbers or use non-uniform grinding, in which the size of the grid cells increases, e.g., logarithmically. I guess it would be problematic if you used the same cell sizes for the ghost and real cell sizes.

I do not think the mentioned approach in partial derivatives can create serious challenges if everything is done truly.

BTW, I guess you can use still a larger dx value :).

I would be happy to answer if there are any other questions :).