# \[Q\]Python and PFC3D give different histories. Why?

**URL:** <https://forum.itascainternational.com/t/q-python-and-pfc3d-give-different-histories-why/602>\
**Category:** PFC\
**Tags:** python\
**Created:** [June 2, 2022, 2:49pm UTC](https://forum.itascainternational.com/t/q-python-and-pfc3d-give-different-histories-why/602 "2022-06-02T14:49:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ppk](https://avatars.discourse-cdn.com/v4/letter/p/7ab992/32.png) [@ppk](https://forum.itascainternational.com/u/ppk)\
**Post date:** [June 2, 2022, 2:49pm UTC](https://forum.itascainternational.com/t/q-python-and-pfc3d-give-different-histories-why/602/1 "2022-06-02T14:49:35Z")

</div>

Simple model which has balls in box.  
Balls are pressed by the wall at the top of the box using servo mechanism.  
The target force exerted on the wall is 100N.

 ![Screen Shot 2022-06-02 at 11.41.20 PM](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/ae1d645566d6d6a2f34f15199aa5defd0f94cc9c.jpeg)

Wall contact force, velocity and displacement are recorded and plotted by Python and PFC3D.

The analysis is performed in deterministic mode.

_ **Question 1 :** _

However, the two methods gives quite different histories for wall velocity and contact force.

This is plot from Python

 ![Screen Shot 2022-06-02 at 10.46.06 PM](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/b015117b15835e636b276b8c51eaa8a62ea5bd73.jpeg)

This is plot from PFC3D

 ![Screen Shot 2022-06-02 at 10.46.25 PM](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/860cc917e4c5af21e0bbc4c14821d8534bf07913.jpeg)

Why the two methods give different results and what is correct one?

_ **Question 2 :** _

According to the sign convention, wall contact force must be positive for this test. However, history of wall contact force obtained from PFC3D command gives negative value. Why?  
The wall velocity history is quite different. Wall velocity history from PFC3D history shows +1.0m/s. It means that the wall is moving upward. This is very confusing results.

![화면 캡처 2022-06-02 161401](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/279d16a9e1b1f0230695c09f9f7d3b87d456b49f.png)

_ **Question 3 :** _

I think that the wall velocity should approach to zero when it reached the target contact force. However, the measured wall velocity still remains maximum velocity limit even when the displacement of the wall reached stationary condition. Is that the wall velocity history is not real measurement?

Results for increased target force (10000N) are as follows.  
Contact force (10500N vs 8200N) and velocity histories are still different.

 ![Screen Shot 2022-06-03 at 12.41.38 AM](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/814f2fb8caed5f8ecff81c253f46dbcaa200dc21.jpeg)

 ![Screen Shot 2022-06-03 at 12.42.23 AM](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/09c7e978f9eb65e39672d5671cc7f1c0dec31746.jpeg)

**"A.py"**

```auto
get_ipython().magic('reset -sf') #command to initialize iPython environment
get_ipython().magic('clear') #iPython screen clear

import itasca as it
import numpy as np
import config as cfg
import importlib
importlib.reload(cfg)

it.command("python-reset-state false") #the model new command will not affect python environment

it.command("""
model new
model configure dynamic
model deterministic off
model large-strain on
model title 'test code'
""")

it.set_threads(96)
print('Program start')
print('No. of threads = ', it.threads())

# SET MODEL DOMAIN
it.set_domain_min(cfg.g_domain_min)
it.set_domain_max(cfg.g_domain_max)

# BALL PROPERTIES
cfg.g_C_MODEL = 'linear' # linear hertz ...
cfg.g_ball_prop_linear = {'density':2700, 'kn':1.0e6, 'ks':1.0e3, 'fric':0.2, 'dp_n':0.7, 'dp_s':0.1}
cfg.g_facet_prop_linear = {'kn':1.0e6, 'ks':1.0e3, 'fric':0.2, 'dp_n':0.7, 'dp_s':0.1}

it.command(f"""
contact cmat default type ball-ball model {cfg.g_C_MODEL} property kn {cfg.g_ball_prop_linear['kn']} \
ks {cfg.g_ball_prop_linear['ks']} fric {cfg.g_ball_prop_linear['fric']} dp_mode 3 dp_nratio {cfg.g_ball_prop_linear['dp_n']} dp_sratio {cfg.g_ball_prop_linear['dp_s']}
contact cmat default type ball-facet model {cfg.g_C_MODEL} property kn {cfg.g_facet_prop_linear['kn']} \
ks {cfg.g_facet_prop_linear['ks']} fric {cfg.g_facet_prop_linear['fric']} dp_nratio {cfg.g_facet_prop_linear['dp_n']} dp_sratio {cfg.g_facet_prop_linear['dp_s']}
""")

# CREATE WALLS
def wall_gen(n, name, p1, p2, p3, p4):
     it.command(f"""
     wall generate id {n} name '{name}' polygon {p1} {p2} {p3} {p4}
     """)

p1, p2, p3, p4 = (-0.1, -0.1, -0.05), (0.1, -0.1, -0.05), (0.1, 0.1, -0.05), (-0.1, 0.1, -0.05)
wall_gen(1, 'bottom', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.1, -0.1, 0.05), (0.1, -0.1, 0.05), (0.1, 0.1, 0.05), (-0.1, 0.1, 0.05)
wall_gen(2, 'top', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, -0.1, -0.05), (-0.08, 0.1, -0.05), (-0.08, 0.1, 0.05), (-0.08, -0.1, 0.05)
wall_gen(3, 'side-1', p1, p2, p3, p4)
p1, p2, p3, p4 = (0.08, -0.1, -0.05), (0.08, 0.1, -0.05), (0.08, 0.1, 0.05), (0.08, -0.1, 0.05)
wall_gen(4, 'side-2', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, -0.08, -0.05), (0.08, -0.08, -0.05), (-0.08, -0.08, 0.05), (0.08, -0.08, 0.05)
wall_gen(5, 'front', p1, p2, p3, p4)
p1, p2, p3, p4 = (-0.08, 0.08, -0.05), (0.08, 0.08, -0.05), (-0.08, 0.08, 0.05), (0.08, 0.08, 0.05)
wall_gen(6, 'back', p1, p2, p3, p4)

# CREATE BALLS
def create_ball(xy_lim, z_lim, r_min, r_max, n_balls):
     it.command(f"""
     ball generate box {-xy_lim} {xy_lim} {-xy_lim} {xy_lim} {-z_lim} {z_lim} number {n_balls} radius {r_min} {r_max} 
     """)
     
create_ball(0.08, 0.05, 0.001, 0.005, 10000)

it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model gravity 9.81
model solve ratio 1e-3
""")

create_ball(0.08, 0.05, 0.001, 0.005, 10000)

it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")

create_ball(0.08, 0.05, 0.001, 0.005, 10000)

it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")

create_ball(0.08, 0.05, 0.001, 0.005, 10000)

it.command(f"""
ball attribute density {cfg.g_ball_prop_linear['density']}
model solve ratio 1e-3
""")

it.add_save_variable('cfg')
it.command(f"""
model save 'gen_balls.sav'
""")

```

**"B.py"**

```auto
get_ipython().magic('reset -sf') #command to initialize iPython environment

import itasca as it
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
it.command("python-reset-state false") #the model new command will not affect python environment
# MODEL CONFIGURATION
it.command("""model new
""")

it.command(f"""
model restore 'gen_balls.sav'
""")

target_force = (0, 0, -100)
top_cap = it.wall.find(2)
top_cap.servo_set_force(target_force)
top_cap.servo_set_vmax(0.1)
top_cap.servo_set_active(True)
print('Target Force (N) is ', top_cap.servo_force())

# HISTORY PARAMETERS
g_ex_1, g_ex_2, g_ex_3, g_ex_4, g_ex_5 = [], [], [], [], []
g_ex_2 = []
g_ex_3 = []
g_ex_4 = []
g_ex_5 = []
#g_ex_6 = []
#g_ex_7 = []

plt.close()

at = np.array([])
at = it.ballarray.pos()

print('Max. z position', np.max(at[:,2]))
max_ball_id = np.argmax(at[:,2]) + 1
print('Max. ball id', max_ball_id) # id
max_ball = it.ball.find(max_ball_id)
print('Max ball position in vector', max_ball.pos())
radius = max_ball.radius()
print('Radius', radius)

target_position = np.max(at[:,2]) + radius * 3
print('Wall target position', target_position)

top_cap.set_pos_z(target_position)

fig = plt.figure(figsize=(15, 6))
gs = gridspec.GridSpec(nrows=2, ncols=2, height_ratios=[1, 1], width_ratios=[1, 1])

def draw_plot(*args):
     global g_ex_1, g_ex_2, g_ex_3
     if it.cycle() % 1000 : return

     g_ex_1.append(it.mech_age())
     g_ex_2.append(top_cap.force_contact_z())
     g_ex_3.append(top_cap.disp_z())
     g_ex_4.append(top_cap.vel_z())
     g_ex_5.append(top_cap.servo_gain())

     plt.figure(1)
     sub_plot_1 = plt.subplot(gs[0])
     plt.xlabel('Time (s)')
     plt.ylabel('Contact Force (N)')
     plt.plot(g_ex_1, g_ex_2, linestyle = "-", color = "Red", linewidth = "0.5")
     plt.grid(True, color = 'Gray', linestyle = '--')
     
     sub_plot_2 = plt.subplot(gs[1])
     plt.xlabel('Time (s)')
     plt.ylabel('Top cap disp. (m)')
     plt.plot(g_ex_1, g_ex_3, linestyle = "-", color = "Black", linewidth = "0.5")

     sub_plot_2 = plt.subplot(gs[2])
     plt.xlabel('Time (s)')
     plt.ylabel('Top cap vel. (m/s)')
     plt.plot(g_ex_1, g_ex_4, linestyle = "-", color = "Black", linewidth = "0.5")
     
     sub_plot_2 = plt.subplot(gs[3])
     plt.xlabel('Time (s)')
     plt.ylabel('Gain')
     plt.plot(g_ex_1, g_ex_5, linestyle = "-", color = "Black", linewidth = "0.5")

     plt.tight_layout()
     plt.draw()
     
it.set_callback("draw_plot",-10.9)

print('Servo gain factor is ', top_cap.servo_gainfactor())
print('Servo gain update is ', top_cap.servo_gainupdate())

it.command(f"""
model deterministic on
history interval 1000
model history name 'd-time' dynamic time-total
wall history name 'z-vel' velocity-z id 2
wall history name 'z-force' force-contact-z id 2
model solve time 0.5
""")

```

_ **"Config.py"** _

```auto
import numpy as np

g_domain_max = (0.1, 0.1, 0.1)
g_domain_min = (-0.1, -0.1, -0.1)
# BALL PROPERTIES
g_C_MODEL = '' # linear hertz ...
g_ball_prop_linear = {}
g_facet_prop_linear = {}

```

---

<div class="post-metadata">

**Author:** ![dblanksma](https://avatars.discourse-cdn.com/v4/letter/d/858c86/32.png) [@dblanksma](https://forum.itascainternational.com/u/dblanksma)\
**Post date:** [June 3, 2022, 6:27pm UTC](https://forum.itascainternational.com/t/q-python-and-pfc3d-give-different-histories-why/602/2 "2022-06-03T18:27:44Z")

</div>

Hello,  
Here are my replies to your questions:  
**Question 1** - In your python function you are getting wall values at cycle location -10.9. Histories will be calculated at the end of the cycle. If you use `it.set_callback("draw_plot",41.0)` (41.0 is right after force-displacement calculations) you will get the same plots between your python plots and PFC plots.

**Question 2** - This is correct. The wall is moving upward and thus creating a tensile (negative) normal force between the balls and the wall. Actually, the wall is fluctuating up and down - see the next comment.

**Question 3** - The wall is actually fluctuating between the two limiting velocities of 0.1 and -0.1. You can see this is you set the history interval to 1. After running the model you provided and changing the history interval to 1, here are the results:

 ![pic](https://canada1.discourse-cdn.com/flex031/uploads/itasca/original/1X/564b3b4df88b52654e987931c5238c4fd92202f9.png)  
If you zoom in you can see each cycle the wall is changing velocity from -0.1 to 0.1 m/s. The same happens for the force. This tells me the servo mechanism probably needs to be changed.

---

<div class="post-metadata">

**Author:** ![ppk](https://avatars.discourse-cdn.com/v4/letter/p/7ab992/32.png) [@ppk](https://forum.itascainternational.com/u/ppk)\
**Post date:** [June 7, 2022, 2:31pm UTC](https://forum.itascainternational.com/t/q-python-and-pfc3d-give-different-histories-why/602/3 "2022-06-07T14:31:05Z")

</div>

Thank you for your kind answer !