Particle Tracking

########################################################################

### This script shows how to use particle tracking APIs in FEFLOW

### A DAC file is loaded, well BC nodes are identified

### and a random well is choosen for demonstration

### Particle tracking is computed using defaults

### Requirement: matplotlib package

########################################################################


import ifm

from ifm import Enum

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.cm import ScalarMappable


fem_file = '..\\femdata\\box_transient.fem'

dac_file = '..\\results\\box_transient.dac'

doc = ifm.loadDocument(dac_file)

pt = doc.createParticleTracer()


# Query all properties of pt:

# print(dir(pt))

# Working with default values:

pt.setDefaults()

# Define a specific porosity:

pt.setPorosityValue(0.25)


# Define settings

pt.setTrackingDirection(1)  # 0 = Forward / 1 = Backward

pt.setTrackingType(0)  # 0 = Only flow velocity / 1 = Random Walk

pt.setTrackingMode(0)  # 0 = steady / 1 = Unsteady


# Generate pathlines based on Well BC coordinates

wellNodes = []

for n in range(0,doc.getNumberOfNodes()):

    if doc.getBcFlowType(n) == Enum.BC_SINGLE_WELL:

        wellNodes.append(n)


# As demonstration we take only one node and its coordinates

wellID = wellNodes[1]

X = doc.getX(wellID)

Y = doc.getY(wellID)

pl = pt.generatePathLine(X, Y)

X_pl = []

Y_pl = []

T_pl = []


# Separate results in lists

for p in pl:

    # print(p)

    X_pl.append(p[0])

    Y_pl.append(p[1])

    T_pl.append(p[2])


### Plotting ###


# Colormap

scales = np.linspace(min(T_pl), max(T_pl), 5)

locs = range(10)

cmap = plt.get_cmap("rainbow") # Reds

norm = plt.Normalize(scales.min(), scales.max())

sm =  ScalarMappable(norm=norm, cmap=cmap)


# Plot pathline

fig, axes = plt.subplots(1,1)

axes.scatter(X_pl, Y_pl, c=T_pl, cmap = cmap, marker='_', linewidth=5, s=120)

axes.scatter(X,Y,c='black',marker='o',linewidth=5, s=120) # plot Well BC node

axes.set_xlim(0,1000)

axes.set_ylim(0,1000)


# Plot colorbar

sm.set_array([])

cbar = fig.colorbar(sm, ax=axes)

cbar.ax.set_title("Time [d]")


# Show plot

plt.show()