Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Question for particle tracking #72

Open
l20000719 opened this issue Jan 8, 2024 · 21 comments
Open

Question for particle tracking #72

l20000719 opened this issue Jan 8, 2024 · 21 comments

Comments

@l20000719
Copy link

l20000719 commented Jan 8, 2024 •
edited
Loading

Dear dfnWorks developers,

Greetings. I have successfully utilized dfnWorks to create a simulation domain with dimensions of 100m x 100m x 100m, generating five sets of fracture networks. Additionally, I've placed a vertical well parallel to the z-axis at the coordinates (0m, 40m, 0m) with a length of 50m and a radius of 4m.

To ensure fractures can cross the well, I've implemented a Python while loop iterating through different seed values. I conducted a fluid flow simulation based on specified boundary conditions in dfn_pressure.in. The pressure is set to 5.001MPa at the y=50m plane and 5MPa at the y=-50m plane.

For particle tracking, I used Python to write data from well_inject.zone into allboundaries.zone, modifying the index to 7. This setup treats the well zone as an inflow boundary in PTDFN_control.dat, with the outflow boundary set as 6 (back_n). Please refer to the attached file (my_DFNcase.zip) for more details.

My question concerns PTDFN_control.dat, which provides 7 options. Although I chose "Option 7. Particles start from in-flow zone nodes positions," I aim for functionality similar to "Option 3." Specifically, I want to release a total of n particles in a specific region (e.g., releasing 3000 particles from the injection well).

I tried "Option 3" and set the region (-4 < x < 4, -36 < y < 44, 0 < z < 50), but result showed no fractures crossed the region and dfnTran wasn't complete, and the code exited. Is there a method to fulfill this requirement, or have I made any mistakes in my configuration?

Best regards,

my_DFNcase.zip
@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

l20000719 commented Jan 8, 2024 •
edited
Loading

Thank you for your prompt reply. The "PTDFN_control.dat" file in the "dfnTrans_input.zip" is my dfnTrans input. It utilizes option 3 to release 3000 particles from the region (-4 < x < 4, -36 < y < 44, 0 < z < 50).
dfnTrans_input.zip

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

Sorry about that, here is the input file as a text file
PTDFN_control.txt

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

l20000719 commented Jan 8, 2024 •
edited
Loading

Alternatively,the dfnTrans input file content is as follows.Sorry for the trouble.

/**********************  INPUT FILES: grid *****************************/
/**** input files with grid of DFN, mainly it's output of DFNGen ******/
param: params.txt
poly: poly_info.dat
inp: full_mesh.inp
stor: full_mesh.stor
boundary: allboundaries.zone
in-flow-boundary: 7
out-flow-boundary: 6

/**************** INPUT FILES: PFLOTRAN flow solution *******************/
PFLOTRAN: yes
PFLOTRAN_vel: darcyvel.dat
PFLOTRAN_cell: cellinfo.dat
PFLOTRAN_uge: full_mesh_vol_area.uge

/**************** INPUT FILES: FEHM flow solution ***********************/
FEHM: no
FEHM_fin: tri_frac.fin

/************************  OUTPUT FILES  ********************************/
out_grid: no
out_3dflow: no
out_init: no
out_tort: no

/*************** output options for particles trajectories ****************/
out_curv: no
out_avs: no
out_traj: no
out_fract: no
out_filetemp: no
allparticles_output: yes

/************* output directories *************************************/
out_dir: dfnTrans_output
out_path: trajectories
out_time: partime.dat

/*************** Intersection Mixing Rule **********************************/
streamline_routing: no

/************************************************************************/
/**************** PARTICLES INITIAL POSITIONS ***************************/
/************************************************************************/

/*Option #1.init_nf: if yes - the same number of particles (init_partn) will be placed 
     on every boundary fracture edge at in-flow boundary, 
     equidistant from each other */
init_nf: no 
init_partn: 10

/*Option #2.init_eqd: if yes - particles will be placed on the same distance from
     each other on all over in-flow boundary edges */  
init_eqd: no  //maximum number of particles that user expects on one fracture edge
init_npart: 10


/*Option #3.all particles start from the same region at in-flow boundary, in a range  
    {in_xmin, in_xmax,in_ymin, in_ymax, in_zmin, in_zmax} */
init_oneregion: yes
in_partn: 3000
in_xmin: -4
in_xmax: 4
in_ymin: 36
in_ymax: 44
in_zmin: 0
in_zmax: 50


/* Option #4.all particles are placed randomly over all fracture surface 
     (not only on boundary edges!) */
init_random: no 
// total number of particles
in_randpart: 50    


/*Option #5.all particles are seed randomly over matrix, 
     they will start travel in DFN from the node/cell that is closest to
     their initial position in rock matrix */     
init_matrix: no
// to obtain these files, run python script RandomPositGener.py
inm_coord: ParticleInitCoordR.dat
inm_nodeID: ClosestNodeR.inp
inm_porosity: 0.02
inm_diffcoeff: 1.0e-12

/* Option #6. particles positions according to in-flow flux weight */
init_fluxw: no
/* Initial number of particles can not be less than 
number of nodes in in-flow boundary. If it is less, 
the number of particles will be increased. */
init_totalnumber: 200 

/* Option #7. Particles start from in-flow zone nodes positions*/
/* for example, nodes defined at well zone/ex file are the position of injection well and fracture intersection */ 
init_well: no
/* Initial number of particles seeded at each node at in-flow zone */ 
init_nodepart: 1

/****************** FLOW AND FRACTURE PARAMETERS **********************/
porosity: 1.0
density: 997.73
satur: 1.0
thickness: 1.0

/************************ APERTURE *********************************/
aperture: yes
aperture_type: frac
aperture_file: aperture.dat

/************* TIME DOMAIN RANDOM WALK ******************************/
tdrw: no
tdrw_porosity: 0.01
tdrw_diffcoeff: 1.0e-9

/********************  TIME ********************************************/
timesteps: 2000000
time_units: years

/**** flux weighted particles*/
flux_weight: yes
seed: 125793

/*********************  Control Plane/Cylinder Output ********************/
ControlPlane: no
control_out: outcontroldir
delta_Control: 10
flowdir: 1

/**************************************************************************/
/endendend/
END

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

In my python code ,I write data from well_inject.zone into allboundaries.zone, modifying the index to 7. Here is a portion of the code.

fzone = open(jobname +"/allboundaries.zone", "r")
lines = fzone.readlines()
a = open(jobname +"/well_inject.zone", "r")
line = a.readlines()
line[1] = '000007     well_inject\n'
writeline = []
for i in lines[:-2]:
    writeline.append(i)
for j in line[1:]:
    writeline.append(j)
with open(jobname +"/allboundaries.zone", "w") as fall:
    fall.writelines(writeline)

This setup treats the well zone as an inflow boundary and set value 7 in PTDFN_control.dat, with the outflow boundary set as 6 (back_n). Here is a full code.

import random
from pydfnworks import *
import os
import numpy as np
example_path = os.getcwd()
jobname =  f"{example_path}/well_example_output"
dfnFlow_file = f"{example_path}/dfn_pressure.in"
dfnTrans_file = f"{example_path}/PTDFN_control.dat"
A = True
run = 1
while A == True:
    DFN = DFNWORKS(jobname,
               dfnFlow_file=dfnFlow_file,
               dfnTrans_file=dfnTrans_file,
               ncpu=20)
    DFN.params['h']['value'] = 0.1
    DFN.params['orientationOption']['value'] = 1
    DFN.params['boundaryFaces']['value'] = [1, 1, 1, 1, 1, 1]
    DFN.params['stopCondition']['value'] = 1
    DFN.params['domainSize']['value'] = [100, 100, 100]
    if run > 1:
        Seed = random.randint(1, 99999)

    else:
        Seed = 3625
    DFN.params['seed']['value'] = Seed

    # Generate the DFN
    rmin_FDMB = 4.5
    rmax_FDMB = 56.4
    FDMB_P32 = [0.15,0.24,0.3,0.1,0.21]
    FDMB_kappa = [20,18,17,16,19]
    FDMB_trend = [65, 344, 281, 174, 175]
    FDMB_plunge = [17, 38, 29, 22, 75]
    powerlaw_a = 1.6
    FDMB_a = 5.3E-11
    FDMB_b = 0.5
    FDMB_P32_total = 0.1 
    for i in range(len(FDMB_P32)):
        DFN.add_fracture_family(shape='rect',
                                distribution='tpl',
                                kappa=FDMB_kappa[i],
                                p32=FDMB_P32_total*FDMB_P32[i],
                                aspect=1.0,
                                trend=FDMB_trend[i],
                                plunge=FDMB_plunge[i],
                                alpha=powerlaw_a,
                                min_radius=rmin_FDMB,
                                max_radius=rmax_FDMB,
                                hy_variable='transmissivity',
                                hy_function='correlated',
                                hy_params={
                                    "alpha": FDMB_a,
                                    "beta": FDMB_b
                                })
    try:
        DFN.make_working_directory(delete=True)
        DFN.check_input()
        DFN.create_network()
    except:
        run += 1
        continue
    # give the region to make sure fracture crossing, if not change the Seed
    xmin, xmax = -10, 10
    ymin, ymax = 30, 50
    zmin, zmax = 0, 50
    data = np.loadtxt(jobname+'/dfnGen_output/intersection_list.dat',
                        skiprows=1, usecols=(2, 3, 4))
    point_in_range = any(xmin <= x <= xmax and
                        ymin <= y <= ymax and
                        zmin <= z <= zmax for x, y, z in data)
    if point_in_range == True:
        # set the well and meshing
        inject_well = {"name": 'inject', "filename": "well_inject.dat","r":4}
        wells = [inject_well]
        os.symlink(f"{example_path}/{inject_well['filename']}",f"{inject_well['filename']}")
        DFN.find_well_intersection_points(wells)
        DFN.mesh_network(well_flag=True)
        DFN.tag_well_in_mesh(wells)
        DFN.cleanup_wells(wells)
        os.chdir(DFN.jobname)
        DFN.combine_well_boundary_zones(wells)
        # combine the well_inject.zone to allboundaries.zone file (inflow zone index number 7)
        fzone = open(jobname +"/allboundaries.zone", "r")
        lines = fzone.readlines()
        a = open(jobname +"/well_inject.zone", "r")
        line = a.readlines()
        line[1] = '000007     well_inject\n'
        writeline = []
        for i in lines[:-2]:
            writeline.append(i)
        for j in line[1:]:
            writeline.append(j)
        with open(jobname +"/allboundaries.zone", "w") as fall:
            fall.writelines(writeline)
        # flow simulation
        DFN.dfn_flow()
        # particle tracking
        DFN.dfn_trans()
        A = False
    else:
        run += 1
        continue
print('Case complete')

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

But when I use the particle tracking option 7 and 1 particles start from in-flow zone for each node , it works! Here is the ouput console

---------------------BOUNDARY CONDITIONS----------------------

 OPEN AND READ FILE: allboundaries.zone 
 
 Read the number and indices of nodes defined in flow in and flow out zones  

 Number of nodes 10903 in flow-in zone.  

 Number of nodes 122 in flow-out  zone.  
Total in-flow volumetric flux =  1.96817e-19 [m^3/s] 
Total out-flow volumetric flux = -1.77831e-11 [m^3/s]
  GRID CHECK starts 
  GRID CHECK is done 
 
 Converting 3d nodes coordinates to 2d xy parallel plane 

----------------VELOCITY RECONSTRUCTION-----------------------

 Darcy's velocities reconstruction 
 Velocities on nodes are calculated 

------------------PARTICLE TRACKING---------------------------

Fracture Intersection Rule: Complete Mixing 
Number of fractures in the in-flow boundary zone: 6



***************************************************
Seeding 10903 particles
There are 2000000 time steps in the current run
***************************************************

***************************************************
Starting Main Loop on Particles
***************************************************
1091 particles out of 10903 have completed (10.01%)
1091 particles out of 1091 have exited successfully.

2181 particles out of 10903 have completed (20.00%)
2181 particles out of 2181 have exited successfully.

3271 particles out of 10903 have completed (30.00%)
3271 particles out of 3271 have exited successfully.

4362 particles out of 10903 have completed (40.01%)
4362 particles out of 4362 have exited successfully.

5452 particles out of 10903 have completed (50.00%)
5452 particles out of 5452 have exited successfully.

6542 particles out of 10903 have completed (60.00%)
6542 particles out of 6542 have exited successfully.

7633 particles out of 10903 have completed (70.01%)
7633 particles out of 7633 have exited successfully.

8723 particles out of 10903 have completed (80.01%)
8723 particles out of 8723 have exited successfully.

9813 particles out of 10903 have completed (90.00%)
9813 particles out of 9813 have exited successfully.

***************************************************
Main Loop on Particles Complete
***************************************************

Number of particles completed 10903, number of particles that went out through out-flow boundary: 5844 


***************************************************
DFNTrans completed successfully.
Finish Time: Mon Jan  8 17:08:48 2024
***************************************************

================================================================================

dfnTrans Complete

Time Required for dfnTrans: 80.66 Seconds

================================================================================

and my allboundary.zone is:

zone
000001     top
nnum
       330
    ...
000002     bottom
nnum
       447
     ...
000003     left_w
nnum
       269
    ...
000004     front_n
nnum
       491
    ...
000005     right_e
nnum
       849
     ...
000006     back_s
nnum
       122
     ...
    396415     396541
000007     well_inject
nnum
     10903
    ...
  
stop

As you can see that the optput console show that the number of nodes in in-flow zone is 10903, and the number of nodes in my rewrited allboudary zone 7 "well_inject" also 10903.

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024

Interesting. Not sure why the box option isn't working then. Can you send me your full_mesh.inp file over transfer.lanl.gov? @jhyman

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

Okay, I have sent it via transfer.lanl.gov, using my email address (l20000719@gmail.com) as the name.

@l20000719
Copy link
Author

image

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024

thanks. I'll look into why option 7 isn't working.

@l20000719
Copy link
Author

l20000719 commented Jan 8, 2024 •
edited
Loading

Apologies for the correction: I intend to release particles in the box region, specifically in the well injection zone, which is option 3, not 7(option 7 works fine).Also thanks for your reply!

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024

ah, you cannot use option 3 on anything but an domain edge. Sorry.

@l20000719
Copy link
Author

l20000719 commented Jan 8, 2024 •
edited
Loading

Okay, so is there any way or other option to let me release the n number of particles from well ?
This is the apprch I can think of： Place the well at the boundary of the domain edge, allowing the release of n particles using option 3 ?

@hymanjd
Copy link
Collaborator

hymanjd commented Jan 8, 2024 via email

@l20000719
Copy link
Author

Thanks a lot!

@l20000719
Copy link
Author

l20000719 commented Jan 9, 2024 •
edited
Loading

I apologize for the interruption again. I have another question. When I mesh the network , I noticed that the mesh around the wells and intersecting fractures does not seem to be refined. (as the figure below).

image

As you can see, the mesh is refined at the intersections of fractures, but there is no refinement between the wells and fractures. Here is the code when I mesh the network.

...
inject_well = {"name": 'inject', "filename": "well_inject.dat","r":4}
wells = [inject_well]
os.symlink(f"{example_path}/{inject_well['filename']}",f"{inject_well['filename']}")
DFN.find_well_intersection_points(wells)
DFN.mesh_network(well_flag=True)
DFN.tag_well_in_mesh(wells)
DFN.cleanup_wells(wells)
os.chdir(DFN.jobname)
DFN.combine_well_boundary_zones(wells)
...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@hymanjd @l20000719