search

vitst / dissolFoam

Solves for flow (Stokes) and transport (steady-state) and moves the mesh according to the reactant flux.
11 stars 2 forks source link

Running the old test cases in dissolFoam-v1912 #5

Open kohoto opened 1 year ago

kohoto commented 1 year ago

Hi, thanks for sharing your amazing work.

I'd like to run dissolFoam with my own geometry and I'm currently studying how to use it. I downloaded OpenFOAM-v1912 and dissoFoam-v1912, and I was able to run the test case.

Now I'd like to run one of the test cases included in the additional information of the 2016 paper (V. Starchenko, C. J. Marra and A. J. C. Ladd, Three-dimensional simulations of fracture dissolution, J. Geophys. Res. Solid Earth, 2016, 121: 6421-6444,) since it is closer to the case I want to run, but I couldn't run it with dissolFoam-v1912.

I'm new to OpenFOAM and I'm not sure which part in the input files I need to change to make it compatible with v1912. I really appreciate it if you could give me brief instructions or references on how the input for 2016 version and v1912 are different.

Thanks a lot!

tonyladd commented 1 year ago

In general the input files are pretty much the same for different versions of OF. Can you post the error messages around the crash. My first guess would be the initial surface - there was code to superpose a random aperture variation on top of the base fracture. We did not use that after about 2018, so it did not get ported to v1912.

Can you post the log

kohoto commented 1 year ago

Hi @tonyladd,

Thanks for your reply. Here are my log files. I'm trying to run case1.

I first ran 'blockMesh' command, and the mesh looked fine to me. After that, I ran ./Allrun command. The program stops after seeing the outputs that are the same as what's in log.reconstruct.

log.decompose.txt log.reconstruct.txt log.run.txt

Thanks for your help!

tonyladd commented 1 year ago

You are going to have to take it step by step through all the stages in the Allrun script. It is crashing on the decompose step which is unusual. I think that may be an out of date library name in your boundary conditions one (or more) of the fields (P, U, C). You should use the names from the v1912 cases. You will have to understand a bit about how OF works to figure out how to merge them. The meshing is an independent step so its not surprising that works You could also try running a serial version - just type dissolFoam. I guess it will crash too, but it if not it will point some odd error with the decompose step.

vitst commented 1 year ago

Hi @kohoto, the problem with the different versions is that over the years we updated the way we code custom boundary conditions.In 2016 we had hardcoded some boundary conditions into a separate libraries which we link to the solver through the controlDict. You can find it in the file: libs ( "libnonLinearFvPatchScalarField.so" "libdanckwertsFvPatchScalarField.so" ); Corresponding files were compiled in libs/danckwerts and libs/nonlinear. In 2019 version we reviewed the practice of creating new boundary condition as a run time compilation coded boundary conditions. That is why the decomposition crashes: it can not find these libraries since they have to be included as coded boundaries, which are compiled automatically during the run. Since you are trying to run fracture dissolution I would recommend to check our updated version which we used in:

V., Starchenko, & Ladd, A. J. C. (2018). The development of wormholes in laboratory-scale fractures: Perspectives from three-dimensional simulations. Water Resources Research, 54, 7946– 7959. https://doi.org/10.1029/2018WR022948

The case can be found in our repository for cases: https://github.com/vitst/dissolCases/tree/master/seededFracture/Pe20lR1

In Zero directory you will find the concentration field file will set 2 coded boundary conditions: for inlet (danckwerts.H) and for walls (linear.H). The corresponding files are located in constant/bcInclude directory. Try to mesh and run this case and let us know if you succeed.

kohoto commented 1 year ago

Hi @tonyladd and @vitst,

Thank you so much for your advice! I understood the reason for crashes. I will try the case you suggested, and I'll read the paper. I will let you know how it goes.

tonyladd commented 1 year ago

OK - the other thing you could try is installing v1706.

kohoto commented 1 year ago

Hi, @vitst and @tonyladd,

I tried the case Pe20lR1, but I still had a problem with loading libnormalMotionSlipPointPatchVectorField.so.

So I compared the input files in cases_1912 and Pe20lR1 and changed the followings, and then the program started to run.

It's still running and I'm not sure what the result is going to look like, but I'd like to ask if these modifications make sense to you. I tried to change the boundary condition from the old version to the new run time compilation-coded boundary conditions.

Thanks for your help!

tonyladd commented 1 year ago

That's good. In general if OF runs the answer is correct (aside from possible convergence issues). It has a lot of error detection built in (which can make it difficult to get running sometimes) but when it does it is very rare that there are errors in the code.

kohoto commented 1 year ago

I see. I hope mine will look the same as yours in the paper.

kohoto commented 1 year ago

Hi @vitst and @tonyladd,

I was able to run the seededFracture case on v1912. But now I have another question I'd like to ask.

I'm confused about the consistency of units.

Why can OF understand the mesh is in 0.1 mm scale and other files in Zero are in normal SI unit?

Also, according to the paper, it looks like the simulation is run until 100 * td (in Eq.(13)), but I'm not sure why OF can recognize the timestep is in the unit of td, not seconds.

Could you tell me how to deal with units?

Thanks in advance!

vitst commented 1 year ago

Hi @kohoto, convertToMeters 1.0 is a OpenFOAM (OF) keyword in blockMeshDict because OF operates in meters. You should look it up in tutorials on dimensional units. However, in our simulation it is not feasible (and not convenient) to operate in units defined by OF. For example, if the unit time is 1 s, the number to represent multiday process is already too large and inconvenient. The dimensionless representation of units is described in our 2016 publication (J. Geophys. Res. Solid Earth, 2016, 121: 6421-6444). You correctly mentioned that unit length now is 0.1 mm or whatever is $h_0$. The description of units you can find in Section 2.2 Dimesionless Equations. Particularly, equation 14 and 16 show length, velocity, and pressure. Equation 18 - concentration. Equation 22 - reaction rate. Finally, equation 24 - time length. In the 2018 paper the concept is pretty much the same. Bottomline: OpenFOAM does not recognize mm or seconds or any other units, but operates with numbers which are defined by the equations coded into solvers. The built-in solvers operate in the units defined by OF developers. Our solver operates in units which we define, which are described in the 2016 paper.

kohoto commented 1 year ago

Hi @vitst,

Thanks for the very clear explanation. Now I understood how the units are dealt with in dissolFoam.

kohoto commented 1 year ago

Hi @vitst and @tonyladd,

I'd like to ask questions about dimensionless numbers again.

  1. I understood how the dimensionless numbers work from your last reply. However, I'm still confused with the characteristic time because $t_d = l_0/k\gamma$ is defined in papers while t0 = l0*A/Q is defined in the comment in transportProperties. Which characteristic time should I use? Or doesn't it matter as long as it's consistent in nu and D?

  2. What is Kinv in transportProperties? At first, I thought it was either k (reaction rate, m/s) or $\kappa$ (dissolution rate constant per unit area of reactive surface, M/m2/s), but the unit is not consistent with Kinv.

Thanks in advance!

vitst commented 1 year ago

Hi @kohoto.

  1. You should follow the the paper definition, ignore comments in transportProperties. We do a lot of experimentation, thus some stuff just slips out.
  2. Kinv is an invert permeability for different set of problems. You should set it to 0 and ignore.

Best

kohoto commented 1 year ago

Hi @vitst,

Thank you so much for your explanation.