Documentation¶

The code is made up of the following modules:

Main Module

Material Points

Material Point Lists

Grid Operators

Constitutive Models

Integrators

Plotting

Main Module¶

The main module. This contains the loop which will be iterated over every timestep.

MPM.main(params)¶

This is the main loop which is repeated every timestep. Currently, this follows the Update Strain Last paradigm (does it?!?).

Parameters: mode (str) – The name of the input file to use.
Returns: int – The return code.

Material Points¶

class particle.Particle(x, y, P, G, p, phi=None)¶

A class defining an individual material point together with functions for mapping to/from the grid

add_to_boundary(P, G)¶

I have no idea what this does.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

energy(P)¶

Calculate energies of this material point for plotting purposes.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

get_basis_functions(P, G)¶

Get the basis functions for this material point.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

get_nodal_forces(P, G)¶

Calculate nearby nodal internal forces and body forces. Tractions are ignored for now.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

get_nodal_mass_momentum(P, G)¶

Pass mass and momentum to nearby nodes.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

get_reference_node(P, G)¶

Find the reference node for this material point.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

increment_gravity(P, G, p)¶

Calculate gravity components for this material point.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

map_gsd_to_mp(P, G)¶

Pass the nodal grainsize information back to the material point.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

map_phi_to_mp(P, G)¶

Pass the nodal grainsize information back to the material point.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

move_material_points(P, G)¶

Update material point position and velocity using nearby nodal mass and momentum.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

move_rigid_points(P, G)¶

Update material point position and velocity for rigid particles with prescribed velocity.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

recover_position(P, G)¶

Find out where this material point should be.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

update_nodal_gsd(P, G)¶

Calculate the grainsize distribution at a node.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

update_strain(P, G)¶

Calculate incremental strain at the particle level from the nodal velocity. Also update the total strain.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

update_stress(P, G, p)¶

Calculate incremental stress at the particle level using the relevant constitutive model.

Parameters

P – A param.Param instance.
G – A grid.Grid instance
p – Which material phase

Material Point Lists¶

class mplist.MatPointList(P, G)¶

A class of material points of a single size. Contains convenience functions for operating on many material points.

cyclic_lr(P, G)¶

Particles which leave the domain in the \(x\) direction are placed back in the other side.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

get_basis_functions(P, G)¶

For every particle, get its basis function.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

get_nodal_forces(P, G)¶

For every particle, get nodal forces. If cyclic boundaries are active, update the nodal cyclic forces.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

get_nodal_gsd(P, G)¶

Measure the volume weighted average size at the node.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

get_nodal_mass_momentum(P, G)¶

For every particle, map mass and momentum to the grid. Then, set momentum at boundaries to zero if applicable.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

get_reference_node(P, G)¶

For every particle, get its reference node.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

inlet_right(P, G)¶

Particles are created at right boundary (\(+x\) direction) of the domain.

This is not physical at all. Needs reimplementation.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

inlet_top(P, G)¶

Particles are created at the top (\(+y\) direction) in the middle of the domain.

This is not physical at all. Needs reimplementation.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

measure_elastic_energy(P, G)¶

I’ve forgotten what this does.

Parameters

P – A param.Param instance.
G – A grid.Grid instance

move_grainsize_on_grid(P, G)¶

Advect grainsize distribution across grid.

Parameters: P – A param.Param instance.

move_material_points(P, G)¶

For every particle, update the particle location.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

outlet_bottom(P, G)¶

Particles are deleted if they move downwards (\(-y\) direction) out of the domain.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

outlet_left(P, G)¶

Particles are deleted if they move left (\(-x\) direction) out of the domain.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

recover_position(P, G)¶

Just used for checking if the code is working. Checks that the particle mapping recovers the real particle position.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

update_forces(P, G)¶

For every particle, update the value of the global gravity.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

update_stress_strain(P, G)¶

For every particle, firstly calculate the incremental strain from nearby nodal velocities. Then use the relevant constitutive law to update the stress at the material point. Finally, if cyclic boundaries are active, update the stress at the boundaries.

Parameters

P – A param.Param instance.
G – A grid.Grid instance.

update_timestep(P)¶

Check through every particle, and if any have yielded, reduce the timestep.

Parameters: P – A param.Param instance.

Grid Operators¶

class grid.Grid(P)¶

This contains all of the methods which operate on the grid directly. The grid is assumed to be a regular lattice, numbered as:

(ny-1)*nx	…	…	…	(ny*nx-1)
…	…	…	…	…
2*nx	…	…	…	(3*nx-1)
nx	…	…	…	(2*nx-1)
0	1	2	…	nx-1

BCs(P)¶: Boundary conditions are applied directly to the external force G.fe

N(x)¶: Calculate shape functions and their derivatives for Lagrange 4-node system

boundary(P)¶

Build the boundary.

This consists of three masks, the same shape as the grid, boundary_v, boundary_h and boundary_tot. These are boundaries that restrain flow in the left/right and up/down directions, respectively. boundary_tot is the sum of the two other boundaries.

This respects the following arguments which may be present in the input file:

has_bottom
has_top
has_left
has_right
no_slip_bottom
wall
box

calculate_gammadot(P, G, smooth=False)¶

Calculate the bulk shear strain rate from the continuum measure of velocity.

Parameters: P – A param.Param instance.

calculate_grad_gammadot(P, G)¶

Calculate the gradient of the absolute value of the shear strain rate using the built-in gradient method.

Parameters: P – A param.Param instance.

calculate_gradient(P, G, Z, smooth=False, verbose=False)¶

Calculate the gradient of any property. Deals with grid points that have no mass (that should’nt contribute to the gradient).

Parameters

P –

A param.Param instance.

Warning

This returns the gradient of the input field … kind of … sometimes

calculate_phi_increment(P)¶

Calculate the incremental change in phi across the grid.

Parameters

P –

A param.Param instance.

Warning

This returns the gradient of the ABSOLUTE VALUE of the shear strain rate

make_cyclic(P, G, params)¶: WORKS (EXCEPT MAYBE WHEN nx < 3) !!!!

update_momentum(P)¶

Update the momentum at every nodal location, as calculated from the nodal change in moment, \(\dot q\). This also…

Parameters: P – A param.Param instance.

update_pk(P, G)¶

Placeholder method until Ebrahim’s model is finished.

Parameters: P – A param.Param instance.

Constitutive Models¶

A set of methods for describing all of the constitutive models.

Each model operates on an individual material point, updating its’ properties. The relevant properties are also mapped to the grid for visualisation.

The models which currently work reliably are:

Rigid
Elastic
Von mises
Newtonian viscosity

constit.bingham(MP, P, G, p)¶

Bingham fluid.

This material model has not been validated but may be functional. YMMV.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.dp(MP, P, G, p)¶

Drucker-Prager \(H^2\) yield criterion.

This material model is ONLY PARTIALLY WORKING - NOT YET VALIDATED!

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.dp_rate(MP, P, G, p)¶

Drucker-Prager yield criterion with rate dependent behaviour.

This material model is PROBALBLY WORKING - NOT YET VALIDATED!

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.elastic(MP, P, G, p)¶

Linear elasticity.

This material model has been validated and is functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.elastic_time(MP, P, G, p)¶

Linear elasticity that is time dependent.

This material model has been validated and is functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.ken_kamrin(MP, P, G, p)¶

\(\mu(I)\) rheology implemented from Ken Kamrin’s MPM paper. Definitely _NOT_ functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.ken_simple(MP, P, G, p)¶

Math: constant mu rheology with stress = 0 when rho < rho_c.

This material model is UNVALIDATED and is DEFINITELY NOT WORKING.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.linear_mu(MP, P, G, p)¶

\(\mu(I)\) rheology implemented only for flowing regime.

This material model has been validated and is PROBABLY functional. (needs more mileage)

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.pouliquen(MP, P, G, p)¶

\(\mu(I)\) rheology implemented only for flowing regime.

This material model has been validated and is PROBABLY functional. (needs more mileage)

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.rigid(MP, P, G, p)¶

Perfectly rigid particles. No stresses are calculated.

This material model has been validated and is functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.viscous(MP, P, G, p)¶

Linear viscosity.

This material model has been validated and is functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

constit.von_mises(MP, P, G, p)¶

Von Mises \(H^2\) yield criterion — see http://dx.doi.org/10.1016/j.ijsolstr.2012.02.003 Eqs 7.11 - 7.15

This material model has been validated and is functional.

Parameters

MP – A particle.Particle instance.
P – A param.Param instance.
G – A grid.Grid instance
p (int) – The particle number.

Integrators¶

A module containing functions for solving partial differential equations on a grid. So far only used for segregation equations.

integrators.div0(a, b)¶

This function replaces nans with zeros when dividing by zero.

Parameters

a – array - Numerator
b – array - Demoniator

Returns

array - a/b, with infs and nans replaced by 0

integrators.limiter(a, b)¶

Pick which limiter to use.

Parameters

a – (array)
b – (array)

Returns

The flux limiting function of a and b

integrators.maxmod(a, b)¶

Maxmod flux limiter.

Parameters

a – (array)
b – (array)

Returns

Maxmod limited function of a and b

integrators.minmod(a, b)¶

Minmod flux limiter.

Parameters

a – (array)
b – (array)

Returns

Minmod limited function of a and b

integrators.superbee(a, b)¶

Superbee flux limiter.

Parameters

a – (array)
b – (array)

Returns

Superbee limited function of a and b

Plotting¶

class plotting.Plotting¶

Class containing all of the plotting functions.

draw_grid(G)¶: Draw an overlay of the grid. Green for both directions, red crosses for horizontal boundaries and red plusses for vertical boundaries.

draw_pretty_grid(G)¶: Draw an overlay of the grid that looks vaguely nice.

savefig(P, name, dpi=150)¶: A method for saving figures in the right place