Merge pull request #114 from ParticulateFlow/release

Release 21.03
This commit is contained in:
Daniel
2021-03-22 16:56:51 +01:00
committed by GitHub
212 changed files with 21477 additions and 1009 deletions

View File

@ -363,6 +363,8 @@ Miscellaneous:
"clear"_clear.html,
"echo"_echo.html,
"extract_surface"_extract_surface.html,
"extrude_surface"_extrude_surface.html,
"if"_if.html,
"include"_include.html,
"jump"_jump.html,
@ -402,8 +404,8 @@ in the command's documentation.
"compute_modify"_compute_modify.html,
"create_atoms"_create_atoms.html,
"create_box"_create_box.html,
"delete_atoms"_delete_atoms.html,
"create_multisphere_clump"_create_multisphere_clump.html,
"delete_atoms"_delete_atoms.html,
"delete_bonds"_delete_bonds.html,
"dielectric"_dielectric.html,
"dihedral_coeff"_dihedral_coeff.html,
@ -413,6 +415,8 @@ in the command's documentation.
"dump"_dump.html,
"dump_modify"_dump_modify.html,
"echo"_echo.html,
"extract_surface"_extract_surface.html,
"extrude_surface"_extrude_surface.html,
"fix"_fix.html,
"fix_modify"_fix_modify.html,
"group"_group.html,
@ -710,6 +714,8 @@ of each style or click on the style itself for a full description:
"ave/atom"_fix_ave_atom.html,
"ave/correlate"_fix_ave_correlate.html,
"ave/euler"_fix_ave_euler.html,
"ave/euler/region"_fix_ave_euler_region.html,
"ave/euler/region/universe"_fix_ave_euler_region.html,
"ave/histo"_fix_ave_histo.html,
"ave/spatial"_fix_ave_spatial.html,
"ave/time"_fix_ave_time.html,
@ -737,6 +743,8 @@ of each style or click on the style itself for a full description:
"evaporate"_fix_evaporate.html,
"execute"_fix_execute.html,
"external"_fix_external.html,
"forcecontrol/region"_fix_forcecontrol_region.html,
"forcecontrol/region/universe"_fix_forcecontrol_region.html,
"freeze"_fix_freeze.html,
"gcmc"_fix_gcmc.html,
"gld"_fix_gld.html,
@ -749,6 +757,8 @@ of each style or click on the style itself for a full description:
"indent"_fix_indent.html,
"insert/pack"_fix_insert_pack.html,
"insert/pack/dense"_fix_insert_pack_dense.html,
"insert/pack/face"_fix_insert_pack_face.html,
"insert/pack/face/universe"_fix_insert_pack_face.html,
"insert/rate/region"_fix_insert_rate_region.html,
"insert/stream"_fix_insert_stream.html,
"insert/stream/moving"_fix_insert_stream_moving.html,
@ -761,6 +771,8 @@ of each style or click on the style itself for a full description:
"lb/viscous"_fix_lb_viscous.html,
"lineforce"_fix_lineforce.html,
"massflow/mesh"_fix_massflow_mesh.html,
"massflow/mesh/face"_fix_massflow_mesh_face.html,
"massflow/mesh/face/universe"_fix_massflow_mesh_face.html,
"meanfreetime"_fix_mean_free_time.html,
"mesh/6dof"_fix_mesh_surface_stress_6dof.html,
"mesh/surface"_fix_mesh_surface.html,

View File

@ -36,6 +36,8 @@ Commands :h1
dump_molfile
dump_custom_vtk
echo
extract_surface
extrude_surface
fix
fix_modify
group

View File

@ -38,7 +38,10 @@ args = list of arguments for a particular style :l
keywords = {output}
{output} values = face or interpolate
dump-identifier = 'stress' or 'id' or 'wear' or 'vel' or 'stresscomponents' or 'owner' or 'area' or 'aedges' or 'acorners' or 'nneigs' or any ID of a "fix mesh/surface"_fix_mesh_surface.html
{euler/vtk} args = none
{euler/vtk} args = zero or more keyword/ value pairs
keywords = {ave_euler} or {cell_center}
{ave_euler} fix-identifier = ID of a fix ave/euler or fix ave/euler/region command
{cell_center} flag = {yes} (default) or {no}
{decomposition/vtk} args = none :pre
{local} args = list of local attributes
@ -382,7 +385,13 @@ commands, you can specify which meshes to dump. If no meshes are specified,
all meshes used in the simulation are dumped.
The {euler/vtk} style dumps the output of a "fix ave/euler"_fix_ave_euler.html
command into a series of VTK files. No further args are expected.
or "fix ave/euler/region"_fix_ave_euler_region.html command into a series of VTK
files. If no further args are specified, the first fix of type ave/euler or
ave/euler/region will be used.
The option {ave_euler} allows to specify the fix-ID of a fix ave/euler or
fix ave/euler/region command.
If {cell_center} is set to {yes}, only the cell centers are dumped, if set to
{no} the cell shapes are dumped.
The {decomposition/vtk} style dumps the processor grid decomposition
into a series of VTK files. No further args are expected.

50
doc/extract_surface.txt Normal file
View File

@ -0,0 +1,50 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
extract_surface command :h3
[Syntax:]
extract_surface inputfile file surfacefile extrusionfile extrude_length l min_rad r mode :pre
inputfile = filename of VTK grid data :ulb,l
file = obligatory keyword :l
surfacefile = filename to write surface VTK data :l
optional extrusion args :l
extrusionfile = filename to write extruded surface VTK data
extrude_length = obligatory keyword
l = total height of extruded volume
min_rad = obligatory keyword
r = distance volume is extruded outwards :pre
mode = output file mode {ascii} or {binary} (optional) :l
:ule
[Examples:]
extract_surface data/grid.vtk file data/surface.vtk data/extrusion.vtk extrude_length 0.1 min_rad 0.025 :pre
[Description:]
Given an unstructured grid VTK file consisting of hexahedral cells,
this command will extract and triangulate the surface and save it to
an unstructured grid VTK file adding a unique ID for each cell face
(i.e. typically two triangles). Furthermore, the surface thus obtained
may be used to create a shell geometry extruding the surface {l-r} distance
units inwards and {r} distance units outwards. If used as insertion volume
for multi-level coarse-grain simulations, {r} should be set to the minimum
particle radius occuring in the simulation. The total extrusion length {l}
should be about 2 particle diameters but is ultimately determined by
the insertion rate and the velocity of the particles.
[Restrictions:]
none
[Default:]
mode = ascii

50
doc/extrude_surface.txt Normal file
View File

@ -0,0 +1,50 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
extrude_surface command :h3
[Syntax:]
extrude_surface inputfile file surfacefile extrusionfile extrude_length l min_rad r mode :pre
inputfile = filename of VTK grid data :ulb,l
file = obligatory keyword :l
surfacefile = filename to write surface VTK data :l
extrusion args :l
extrusionfile = filename to write extruded surface VTK data
extrude_length = obligatory keyword
l = total height of extruded volume
min_rad = obligatory keyword
r = distance volume is extruded outwards :pre
mode = output file mode {ascii} or {binary} (optional) :l
:ule
[Examples:]
extrude_surface data/quadsurface.vtu file data/trisurface.vtu data/extrusion.vtk extrude_length 0.1 min_rad 0.025 :pre
[Description:]
Given an unstructured grid dataset or polydata surface VTK file, this command
will triangulate the surface cells (polygon, pixel or quad) and save it to an
unstructured grid VTK file adding a unique ID for each cell face (i.e. typically
two triangles).
Furthermore, the surface thus obtained is used to create a 3D geometry extruding
the surface {l-r} distance units inwards and {r} distance units outwards.
If used as insertion volume for multi-level coarse-grain simulations,
{r} should be set to the minimum particle radius occuring in the simulation.
The total extrusion length {l} should be about 2 particle diameters but is
ultimately determined by the insertion rate and the velocity of the particles.
[Restrictions:]
none
[Default:]
mode = ascii

View File

@ -0,0 +1,115 @@
"LIGGGHTS WWW Site"_liws - "LAMMPS WWW Site"_lws - "LIGGGHTS Documentation"_ld - "LIGGGHTS Commands"_lc :c
:link(liws,http://www.cfdem.com)
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix ave/euler/region command :h3
fix ave/euler/region/universe command :h3
[Syntax:]
fix ID group-ID ave/euler/region nevery N region reg-ID keywords values
fix ID group-ID ave/euler/region/universe nevery N region reg-ID keywords values ukeywords uvalues :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
ave/euler/region = style name of this fix command :l
nevery = obligatory keyword :l
N = calculate average values every this many timesteps (also sending interval in universe version) :l
region = obligatory keyword :l
reg-ID = ID of region with style mesh/hex defining the grid :l
zero or more keyword/value pairs may be appended :l
keyword = {basevolume_region} :l
{basevolume_region} value = region-ID
region-ID = correct grid cell volume based on this region :pre
one or more ukeyword/uvalue pairs must be appended for the universe version of this command :l
ukeywords = {send_to_partition} (obligatory) or {sync} :l
{send_to_partition} value = partition
partition = partition to send data to in multi-partition simulations
{sync} value = mode
mode = {yes} to use MPI_Ssend, {no} to use MPI_Bsend for communication between partitions :pre
:ule
[Examples:]
fix 1 all ave/euler/region nevery 10 region ave_reg
fix 2 all ave/euler/region/universe nevery 10 region ave_reg send_to_partition 2 sync yes :pre
[Description:]
Calculate cell-based averages of velocity, radius, volume fraction,
and pressure (-1/3 * trace of the stress tensor) every few timesteps,
as specified by the {nevery} keyword. The cells are taken from the
specified region.
Note that velocity is favre (mass) averaged, whereas radius is arithmetically
averaged. To calculate the stress, this command internally uses a
"compute stress/atom"_compute_stress_atom.html . It includes the convective
term correctly for granular particles with non-zero average velocity
(which is not included in "compute stress/atom"_compute_stress_atom.html).
However, it does not include bond, angle, diahedral or kspace contributions
so that the stress tensor finally reads
:c,image(Eqs/stress_tensor_granular.png)
where vave is the (cell-based) average velocity.
The first term is a kinetic energy contribution for atom {I}. The
second term is a pairwise energy contribution where {n} loops over the
{Np} neighbors of atom {I}, {r1} and {r2} are the positions of the 2
atoms in the pairwise interaction, and {F1} and {F2} are the forces on
the 2 atoms resulting from the pairwise interaction.
The {basevolume_region} option allows to specify a region that
represents the volume which can theoretically be filled with
particles. This will then be used to correct the basis of the averaging
volume for each cell in the grid. For example, if you use a cylindrical
wall, it makes sense to use an identical cylindrical region for
the {basevolume_region} option, and the command will correctly
calculate the volume fraction in the near-wall cells.
the calculation of overlap between grid cells and the region
is done using a Monte-Carlo approach.
If LIGGGHTS is invoked in multi-partition mode (cf. "Command-line options"_Section_start.html#start_7),
where each partition represents a separate coarse-graining level of the multi-level
coarse-grain model, the universe version of this command needs to be used to
provide input data for a fix of style {forcecontrol/region/universe} on the
coupled partition.
The coupled partition must be specified via the {send_to_partition} option.
Furthermore, in this mode the {nevery} option also specifies the sending interval.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart
files"_restart.html. None of the "fix_modify"_fix_modify.html options
are relevant to this fix.
This fix computes the above-mentioned quantities for output via a
"dump euler/vtk"_dump.html command. The values can
only be accessed on timesteps that are multiples of {nevery} since that
is when calculations are performed.
No parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command. This fix is not invoked during "energy
minimization"_minimize.html.
[Restrictions:]
Volume fractions and stresses are calculated based on the specified
grid, so volume fractions and stresses near walls that are not
aligned with the grid will be incorrect.
[Related commands:]
"compute"_compute.html, "compute stress/atom"_compute_stress_atom.html,
"fix ave/atom"_fix_ave_atom.html, "fix ave/euler"_fix_ave_euler.html,
"fix ave/histo"_fix_ave_histo.html, "fix ave/time"_fix_ave_time.html,
"fix ave/spatial"_fix_ave_spatial.html, "partition"_partition.html
[Default:] none

View File

@ -0,0 +1,145 @@
"LIGGGHTS WWW Site"_liws - "LAMMPS WWW Site"_lws - "LIGGGHTS Documentation"_ld - "LIGGGHTS Commands"_lc :c
:link(liws,http://www.cfdem.com)
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix forcecontrol/region command :h3
fix forcecontrol/region/universe command :h3
[Syntax:]
fix ID group-ID forcecontrol/region control_keywords control_values
fix ID group-ID forcecontrol/region/universe control_keywords control_values ukeywords uvalues :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
forcecontrol/region = style name of this fix command :l
control_keyword/value pairs :l
control_keywords = {ctrlPV} (obligatory) or {actual_val} (obligatory) or {target_val} (obligatory) or {velocity_limit} or {kp} or {ki} or {kd} or {actual_cg} or {target_cg} or {max_cg} :l
{ctrlPV} values = {stress} or {velocity}
stress = use stress as control process value
velocity = use velocity as control process value
{actual_val} values = fix-ID
fix-ID = ID of a "fix ave/euler/region"_fix_ave_euler_region.html command that holds the actual values for the controller
{target_val} values = fix-ID
fix-ID = ID of a "fix ave/euler/region"_fix_ave_euler_region.html or "fix ave/euler/region/universe"_fix_ave_euler_region.html command that holds the target values for the controller
{velocity_limit} values = mode value
mode = {off} or {on} or {relative} or {absolute}
value = N/A to {off} or {on} mode, allowed deviation from target velocity for {relative} and {absolute} mode
{kp} values = k
k = proportional constant for PID controller
{ki} values = k
k = integral constant for PID controller
{kd} values = k
k = differential constant for PID controller
{actual_cg} values = cg
cg = coarse-grain ratio of the actual particles (N/A to universe version)
{target_cg} values = cg
cg = coarse-grain ratio of the target particles (N/A to universe version)
{max_cg} values = cg
cg = maximum coarse-grain ratio in the simulation (N/A to universe version) :pre
two ukeyword/uvalue pairs must be appended for the universe version of this command :l
ukeywords = {receive_from_partition} and {couple_every} :l
{receive_from_partition} value = partition
partition = partition to receive data from in multi-partition simulations
{couple_every} value = interval
interval = time interval for receiving data :pre
:ule
[Examples:]
fix stressctrl resolved forcecontrol/region ctrlPV stress actual_val resolvedave target_val coarseave kp 1. ki 10000. cg 2
fix velocityctrl coarse forcecontrol/region ctrlPV velocity actual_val coarseave target_val resolvedave kp 0.1 :pre
[Description:]
This fix applies cell-wise forces on particles inside a hexhedral grid
until either a granular pressure (for {ctrlPV} = stress) or an average
velocity (for {ctrlPV} = velocity) is reached.
The actual cell-wise stress or velocity values are calculated by a
"fix ave/euler/region"_fix_ave_euler_region.html, which is defined via
the keyword {actual_val}.
The cell-wise target stress or velocity values are calculated by a
"fix ave/euler/region"_fix_ave_euler_region.html, which is defined via
the keyword {target_val}.
For {ctrlPV} = stress the controller force is applied only in a boundary
layer of a cell with a thickness of 1.3 times the maximum particle diameter
d_max. After this layer the force decays in a sine squared fashion within
0.25*d_max. The direction of the force is determined via the grid geometry
which should hold a vectorial element property with the label 'stress_ctrl_dir'
(e.g. by specifying a cell_data property in the VTK geometry file).
For {ctrlPV} = velocity the controller force is kept constant over
the whole cell.
The controller itself is a proportional-integral-derivative (PID)
controller which is defined by 3 constants {kp}, {ki}, {kd}:
output = kp * error + ki * errorsum + kd * errorchange
where 'error' is the current deviation of the control process value to the target value,
'errorsum' is the time integration (sum) of the error values and 'errorchange' its derivative.
In difficult to control systems (i.e. hard to reach stress target values) the
controller force may produce undesirably high particle velocities.
In this case the {velocity_limit} option may be used to bring velocities of the
controlled particles (or more precisely the controller force) in line with the
velocity range of the target system.
The option {velocity_limit on} will try to limit particle velocities to the
target system's max/min velocities times (1 plus/minus 0.2), i.e. plus/minus 20%.
Similarly, the option {velocity_limit relative value} will try to limit
particle velocities to the target system's max/min velocities times (1 plus/minus {value}).
The option {velocity_limit absolute value} will try to limit particle velocities
to the target system's max/min velocities plus/minus {value}.
The {actual_cg} parameter allows to specify the coarse-grain ratio of the
controlled group of particles in single-partition multi-level coarse-grain
simulations.
In multi-partition simulations this coarse-grain ratio is automatically deduced
from the "coarsegraining"_coarsegraining.html command active in the same
partition.
The {target_cg} parameter allows to specify the coarse-grain ratio of the target
particles in single-partition multi-level coarse-grain simulations.
In multi-partition simulations this coarse-grain ratio is automatically deduced
from the "coarsegraining"_coarsegraining.html command active in the coupled
partition.
The {max_cg} parameter allows to specify the maximum coarse-grain ratio of
particles in single-partition multi-level coarse-grain simulations.
This coarse-grain ratio is need to properly scale the controlled region of
stress-based controllers.
If LIGGGHTS is invoked in multi-partition mode (cf. "Command-line options"_Section_start.html#start_7),
where each partition represents a separate coarse-graining level of the multi-level
coarse-grain model, the universe version of this command needs to be used.
In this case the coupled partition must be specified via the {receive_from_partition}
option.
Furthermore, {target_val} needs to specify a fix of style {ave/euler/region/universe}
on the coupled partition and the interval for receiving data from that fix needs
to be set via the {couple_every} option.
[Restart, fix_modify, output, run start/stop, minimize info:]
This fix supports "fix_modify"_fix_modify.html with option {activate}/{deactivate} = cell-ID
to turn on/off the force controller for individual grid cells.
This fix also supports "fix_modify"_fix_modify.html with option {massflow_correction_on}/{massflow_correction_off} = cell-ID
to turn on/off corrections to a 'velocity' controller to achieve proper mass flow rates.
[Related commands:]
"fix ave/euler/region"_fix_ave_euler_region.html, "region mesh/hex"_region.html,
"partition"_partition.html
[Default:]
{velocity_limit} = {off}, {kp} = 1, {ki} = 0, {kd} = 0, {target_cg} = 1,
{actual_cg} = taken from {coarsegraining} command or 1, {max_cg} = {actual_cg}

View File

@ -0,0 +1,163 @@
"LIGGGHTS WWW Site"_liws - "LAMMPS WWW Site"_lws - "LIGGGHTS Documentation"_ld - "LIGGGHTS Commands"_lc :c
:link(liws,http://www.cfdem.com)
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix insert/pack/face command :h3
fix insert/pack/face/universe command :h3
[Syntax:]
fix ID group-ID insert/pack/face seed seed_value general_keywords general_values pack_face_keywords pack_face_values
fix ID group-ID insert/pack/face/universe seed seed_value general_keywords general_values pack_face_keywords pack_face_values ukeyword uvalue :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
insert/pack/face = style names of this fix command :l
seed = obligatory keyword :l
seed_value = random # seed (positive integer) :l
one or more general keyword/value pairs can be appended :l
general_keywords = {verbose} or {maxattampt} or {insert_every} or {all_in} :l
{verbose} = yes or no
{maxattempt} value = ma
ma = max # of insertion attempts per atom (positive integer)
{insert_every} value = ie
ie = every how many time-steps particles are inserted (also receiving interval in universe version) - insertion happens periodically (positive integer)
{start} value = ts
ts = time-step at which insertion should start (positive integer larger than current time-step)
{all_in} value = yes or no :pre
following the general keyword/value section, one or more pack_face keyword/value pairs can be appended for the fix insert/pack/face command :l
pack_face_keywords = {region} or {massflow_face} or {cg} or {type_offset} or {ntry_mc} :l
{region} value = region-ID
region-ID = ID of the region mesh/hex where the particles will be generated
{massflow_face} values = fix-ID
fix-ID = ID of a "fix massflow/mesh/face"_fix_massflow_mesh_face.html or "fix massflow/mesh/face/universe"_fix_massflow_mesh_face.html command
{cg} value = cg
cg = coarse grain ratio of particles to insert (N/A to universe version)
{type_offset} value = to
to = offset in atom type compared to particles measured by "fix massflow/mesh/face"_fix_massflow_mesh_face.html
{ntry_mc} values = n
n = number of Monte-Carlo steps for calculating the region's volume (positive integer)
{temperature} value = {yes} or {no}
yes = to set the temperature of the inserted particles
{chemistry} value = {yes} or {no}
yes = to set the reduction state of the inserted particles :pre
one ukeyword/uvalue pair must be appended for the universe version of this command :l
ukeywords = {receive_from_partition} :l
{receive_from_partition} value = partition
partition = partition to receive data from in multi-partition simulations :pre
:ule
[Examples:]
fix ins resolved insert/pack/face seed 1001 maxattempt 500 insert_every 500 all_in yes region hexregion massflow_face massflow cg 1 type_offset 1 ntry_mc 10000 :pre
[Description:]
Insert particles recorded by "fix massflow/mesh/face"_fix_massflow_mesh_face.html
into a granular run every few timesteps within the specified region,
as defined via the {region} keyword. The region of type mesh/hex must
hold the element property 'face_id' corresponding to the 'face_id' property
of the mesh used in fix massflow/mesh/face. This command uses a distributiontemplate
generated automatically from the data collected by fix massflow/mesh/face to define
the properties of the inserted particles.
The {verbose} keyword controls whether statistics about particle
insertion is output to the screen each time particles are inserted.
At each insertion step, fix insert/pack/face tries to insert all particles
recorded by fix massflow/mesh/face since the last insertion.
The frequency of the particle insertion can be controlled by the
keyword {insert_every}, which defines the number of time-steps between
two insertions.
The {start} keyword can be used to set the time-step at which the insertion
should start.
Inserted particles are assigned the atom types recorded by fix massflow/mesh/face
which can be altered via the {type_offset} keyword. The particles are assigned to
2 groups: the default group "all" and the group specified in the fix insert command.
Overlap is checked for at insertion, both within the inserted particle package
and with other existig particles. The number of insertion attempts per particle
can be specified via the {maxattempt} keyword. Each timestep particles are inserted,
the command will make up to a total of M tries to insert the new particles without
overlaps, where M = # of inserted particles * {maxattempt}.
If unsuccessful at completing all insertions, a warning will be printed.
The {all_in} flag determines if the particle is completely contained
in the insertion region ({all_in yes}) or only the particle center
({all_in no}).
The initial velocity is determined by the face- and mass-averaged data
gathered by fix massflow/mesh/face.
[Description for fix insert/pack/face:]
This command must use the {region} keyword to define an insertion
volume. The specified region must have been previously defined with a
"region mesh/hex"_region.html command. Dynamic regions are not supported
as insertion region. Each timestep particles are inserted, they are placed
randomly inside the hexahedral cell of the insertion volume that corresponds
to the face of the surface mesh used by fix massflow/mesh/face.
The {massflow_face} must be used to specify a
"fix massflow/mesh/face"_fix_massflow_mesh_face.html command that
determines the number and parameters of the particles to insert.
The {cg} value determines the coarse grain ratio of the particles to insert.
The {type_offset} option can be used to change the atom type of the
inserted particles from the atom type of the particles recorded by
fix massflow/mesh/face. The given offset is added to the original atom type.
To determine the volume of each hexahedron of the insertion region,
a Monte Carlo approach is used. The {ntry_mc} keyword is used to control
the number of MC tries that are used for the volume calculation.
If LIGGGHTS is invoked in multi-partition mode (cf. "Command-line options"_Section_start.html#start_7),
where each partition represents a separate coarse-graining level of the multi-level
coarse-grain model, the universe version of this command needs to be used.
In this case the coupled partition must be specified via the {receive_from_partition}
option.
Furthermore, {massflow_face} needs to specify a fix of style {massflow/mesh/face/universe}
on the coupled partition.
In this mode the {nevery} option also specifies the interval for receiving data.
[Restart, fix_modify, output, run start/stop, minimize info:]
Information about this fix is written to "binary restart
files"_restart.html. This means you can restart a simulation
while inserting particles, when the restart file was written during the
insertion operation.
None of the "fix_modify"_fix_modify.html options are relevant to this
fix. A global vector is stored by this fix for access by various "output
commands"_Section_howto.html#howto_15. The first component of the vector is the
number of particles already inserted, the second component is the mass
of particles already inserted. No parameter of this fix can be
used with the {start/stop} keywords of the "run"_run.html command.
This fix is not invoked during "energy minimization"_minimize.html.
[Restrictions:]
Dynamic regions are not supported as insertion region.
[Related commands:]
"fix insert/pack"_fix_insert_pack.html, "fix insert/rate/region"_fix_insert_rate_region.html,
"fix insert/stream"_fix_insert_stream.html, "fix deposit"_fix_deposit.html,
"fix pour"_fix_pour.html", region"_region.html, "partition"_partition.html
[Default:]
The defaults are {maxattempt} = 50, {all_in} = {no}, {start} = next time-step,
{cg} = 1, {type_offset} = 0, {ntry_mc} = 100000

View File

@ -0,0 +1,146 @@
"LIGGGHTS WWW Site"_liws - "LAMMPS WWW Site"_lws - "LIGGGHTS Documentation"_ld - "LIGGGHTS Commands"_lc :c
:link(liws,http://www.cfdem.com)
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
fix massflow/mesh/face command :h3
fix massflow/mesh/face/universe command :h3
[Syntax:]
fix ID group-ID massflow/mesh mesh mesh-ID vec_side vx vy vz keywords values
fix ID group-ID massflow/mesh/universe mesh mesh-ID vec_side vx vy vz keywords values ukeywords uvalues :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
massflow/mesh/face = style name of this fix command :l
mesh = obligatory keyword :l
mesh-ID = ID of a "fix mesh/surface"_fix_mesh_surface.html command :l
zero or more keyword/value pairs may be appended to args :l
keywords = {check_property} or {count} or {point_at_outlet} or {append} or {file} or {screen} or {delete_atoms} :l
{check_property} property value
property = name of atom property to check
value = reference atom property value to check against
{count} value = {once} or {multiple}
once = count particles only once
multiple = allow particles to be counted multiple times
{point_at_outlet} pointX pointY pointZ
pointX pointY pointZ = coordinates of point on the outlet side of the surface
{inside_out}
use this in connection with {point_at_outlet} to flip direction particle counting
{file} value = filename
{append} value = filename
filename = name of the file to print diameter, position and velocity values of the particles
{screen} value = {yes} or {no}
{delete_atoms} value = {yes} or {no}
yes = to remove the particles that pass through the mesh surface
{cg} value = cg
cg = coarse grain ratio of measured particles (N/A to universe version)
{temperature} value = {yes} or {no}
yes = to measure the temperature of the particles
{chemistry} value = {yes} or {no}
yes = to measure the reduction state of the particles :pre
two ukeyword/uvalue pairs must be appended for the universe version of this command :l
ukeywords = {send_to_partition} and {couple_every} :l
{send_to_partition} value = partition
partition = partition to send data to in multi-partition simulations
{couple_every} value = interval
interval = time interval for sending data :pre
:ule
[Examples:]
fix massflowcg2t all massflow/mesh/face mesh surfacet count once check_property i_cglevel 2 cg 2 inside_out :pre
fix massflowcg2b all massflow/mesh/face mesh surfaceb count multiple check_property i_cglevel 2 cg 2 file post/massflowprop.txt :pre
[LIGGGHTS vs. LAMMPS Info:]
This LIGGGHTS command is not available in LAMMPS.
[Description:]
Fix massflow/mesh/face tracks how many particles penetrate through each face of
a mesh surface, as defined by a "fix mesh/surface"_fix_mesh_surface.html command.
It counts the total number of particles and the associated mass. It calculates the
averaged velocity of the crossing particles and stores their size distribution as
well as their atom types. Only particles part of {group} are eligible for counting.
If {check_property} is specified, only particles with the given property value are
considered. {check_property} may refer to a "fix poperty/atom"_fix_property.html
or a fix property/atom/lammps.
Particles are counted if they cross from the inner side of the mesh to the outer
side of the mesh. The outer side can be defined by the winding order of the mesh
vertices or by specifying a point at the outlet side of the mesh (keyword {point_at_outlet}).
The following restriction applies in case {point_at_outlet} is used: the {count}
value has to be set to once.
The keyword {point_at_outlet} is especially useful in case a cylindrically-shaped
surface is used. The {point_at_outlet} value should be on the cylinder axis in
this case. If you like to track particles moving away from the cylinder axis,
specify the {point_at_outlet} on the axis, and use the keyword
{inside_out} to flip the direction.
If {count} = once, then each particle is only counted once, for {count} = multiple
a particle contributes to the statistics each time it crosses the mesh face.
This can happen e.g. in the case of periodic boundary conditions or in re-circulating
flow conditions.
The timestep of the crossing, the ID of the crossed face, the atom ID, as well as
the diameter, position and velocity of the particles can be written into a file
using the {file} keyword and specifying a filename.
If the {screen} keyword is used, output by this fix to the screen and
logfile can be turned on or off as desired.
If the {delete_atoms} keyword is used then the particles passing through the mesh
surface are deleted at the next re-neighboring step.
If LIGGGHTS is invoked in multi-partition mode (cf. "Command-line options"_Section_start.html#start_7),
where each partition represents a separate coarse-graining level of the multi-level
coarse-grain model, the universe version of this command needs to be used to
provide input data for a fix of style {insert/pack/face/universe} on the
coupled partition.
The coupled partition must be specified via the {send_to_partition} option.
Furthermore, in this mode the {couple_every} option specifies the sending interval.
[Restart, fix_modify, output, run start/stop, minimize info:]
Information about this fix is written to "binary restart files"_restart.html .
This fix computes a per-atom vector (the marker) which can be accessed by
various "output commands"_Section_howto.html#howto_15. The per-atom vector
(i.e., the marker) can be accessed by dumps by using "f_massflowface_ID".
This fix also computes a global vector of length 10. This vector can be
accessed via "f_ID", where ID is the fix id. The first vector component
is equal to the total mass which has crossed the mesh surface, the second vector
component indicates the (resolved) particle count. The third vector component
is equal to the total mass which has crossed the mesh surface since the last output
divided by the time since the last output (i.e., the mass flow rate), the fourth vector
component indicates the (resolved) particle count since the last output divided by the time
since the last output (i.e., the number rate of particles). The fifth and sixth vector
components are the deleted mass and the number of deleted particles. The seventh component
vector component indicates the (resolved) particle count since the last output.
The last three vector components are the favre-averaged velocity components.
This vector can also be accessed by various "output commands"_Section_howto.html#howto_15.
[Restrictions:]
The mesh surface must provide the element property 'face_id', e.g. by loading it from
a VTK file using the {cell_data} option of the "fix mesh/surface"_fix_mesh_surface.html
command. All triangles with the same ID are considered to form a face.
[Related commands:]
"compute nparticles/tracer/region"_compute_nparticles_tracer_region.html,
"fix massflow/mesh"_fix_massflow_mesh.html, "partition"_partition.html
[Default:]
{count} = multiple, {inside_out} = false, {delete_atoms} = false, {cg} = 1

View File

@ -16,18 +16,19 @@ fix mesh/surface/omp command :h3
fix ID group-ID mesh/surface file filename premesh_keywords premesh_values mesh_keywords mesh_values surface_keywords surface_values
fix ID group-ID mesh/surface/planar file filename premesh_keywords premesh_values mesh_keywords mesh_values surface_keywords surface_values :pre
ID, is documented in "fix"_fix.html command. :ulb,l
ID is documented in "fix"_fix.html command. :ulb,l
mesh/surface or mesh/surface/planar = style name of this fix command :l
file = obligatory keyword :l
filename = name of STL or VTK file containing the triangle mesh data :l
zero or more premesh_keywords/premesh_value pairs may be appended :l
premesh_keyword = {type} or {precision} or {heal} or {element_exclusion_list} or {verbose} :l
premesh_keyword = {type} or {precision} or {heal} or {element_exclusion_list} or {cell_data} or {verbose} :l
{type} value = atom type (material type) of the wall imported from the STL file
{precision} value = length mesh nodes this far away at maximum will be recognized as identical (length units)
{heal} value = auto_remove_duplicates or no
{element_exclusion_list} values = mode element_exlusion_file
mode = read or write
element_exlusion_file = name of file containing the elements to be excluded
{cell_data} value = yes or no to load per element data from VTK files
{verbose} value = yes or no :pre
zero or more mesh_keywords/mesh_value pairs may be appended :l
mesh_keyword = {scale} or {move} or {rotate} or {temperature} :l

View File

@ -16,9 +16,10 @@ fix ID group-ID nve/sphere :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
nve/sphere = style name of this fix command :l
zero or more keyword/value pairs may be appended :l
keyword = {update} :l
keyword = {update} or {implicit_integration} :l
{update} value = {dipole}
dipole = update orientation of dipole moment during integration :pre
dipole = update orientation of dipole moment during integration
{implicit_integration} value = {yes} or {no} :pre
:ule
[Examples:]
@ -42,6 +43,10 @@ during the time integration. This option should be used for models
where a dipole moment is assigned to particles via use of the
"atom_style dipole"_atom_style.html command.
{implicit_integration} (default {no}) treats the drag force in a coupled simulation implicitly
if the same option is activated in fix/cfd/coupling/force/implicit and the CFDEM
drag force model.
:line
Styles with a {cuda}, {gpu}, {omp}, or {opt} suffix are functionally

View File

@ -17,11 +17,12 @@ print = style name of this fix command :l
N = print every N steps :l
string = text string to print with optional variable names :l
zero or more keyword/value pairs may be appended :l
keyword = {file} or {append} or {screen} or {title} :l
keyword = {file} or {append} or {screen} or {title} or {starttime} or {endtime} :l
{file} value = filename
{append} value = filename
{screen} value = {yes} or {no}
{title} value = string
{start/endtime} value = float
string = text to print as 1st line of output file :pre
:ule
@ -65,7 +66,10 @@ printed as the first line of the output file, assuming the {file} or
# Fix print output for fix ID :pre
where ID is replaced with the fix-ID. If the specified string equals
{none}, no title line is printed
{none}, no title line is printed.
{starttime} and/or {endtime} allow to specifiy a time window when output
is printed.
[Restart, fix_modify, output, run start/stop, minimize info:]

View File

@ -13,6 +13,7 @@ Fixes :h1
fix_ave_atom
fix_ave_correlate
fix_ave_euler
fix_ave_euler_region
fix_ave_histo
fix_ave_spatial
fix_ave_time
@ -40,6 +41,7 @@ Fixes :h1
fix_evaporate
fix_execute
fix_external
fix_forcecontrol_region
fix_freeze
fix_gcmc
fix_gld
@ -51,6 +53,7 @@ Fixes :h1
fix_indent
fix_insert_pack
fix_insert_pack_dense
fix_insert_pack_face
fix_insert_rate_region
fix_insert_stream
fix_insert_stream_moving
@ -64,6 +67,7 @@ Fixes :h1
fix_limit_vel
fix_lineforce
fix_massflow_mesh
fix_massflow_mesh_face
fix_mean_free_time
fix_mesh_surface
fix_mesh_surface_stress

View File

@ -37,7 +37,7 @@ The coefficient of rolling friction (rmu) must be defined as
fix id all property/global coefficientRollingFriction peratomtypepair n_atomtypes value_11 value_12 .. value_21 value_22 .. .
(value_ij=value for the coefficient of rolling friction between atom type i and j; n_atomtypes is the number of atom types you want to use in your simulation) :pre
This coefficient rmu is equal to the rmu as defined in the "CDT model"_gran_rolling_friction_cdt.html. :pre
This coefficient rmu is equal to the rmu as defined in the "CDT model"_gran_rolling_friction_cdt.html.
IMPORTANT NOTE: You have to use atom styles beginning from 1, e.g. 1,2,3,...

View File

@ -27,7 +27,7 @@ partition yes 6* fix all nvt temp 1.0 1.0 0.1 :pre
This command invokes the specified command on a subset of the
partitions of processors you have defined via the -partition
command-line switch. See "Section_start 6"_Section_start.html#start_7
command-line switch. See "Getting Started"_Section_start.html#start_7
for an explanation of the switch.
Normally, every input script command in your script is invoked by

View File

@ -39,6 +39,17 @@ style = {delete} or {block} or {cone} or {cylinder} or {plane} or {prism} or {sp
{sphere} args = x y z radius
x,y,z = center of sphere (distance units)
radius = radius of sphere (distance units)
{mesh/hex} args = file filename scale s move offx offy offz rotate phix phiy phiz cell_data cd
file = obligatory keyword
filename = name of ASCII VTK file containing the VTK hex-mesh data
scale = obligatory keyword
s = scaling factor for the mesh (dimensionless)
move = obligatory keyword
offx,offy,offz = offset for the mesh (distance units)
rotate = obligatory keyword
phix,phiy,phiz = angle of mesh rotation around x-, y-, and z-axis (in grad)
cell_data = obligatory keyword
cd = {yes} or {no} to load cell properties from VTK file
{mesh/tet} args = file filename scale s move offx offy offz rotate phix phiy phiz
file = obligatory keyword
filename = name of ASCII VTK file containing the VTK tet-mesh data
@ -217,6 +228,12 @@ for a geometric description of triclinic boxes, as defined by LAMMPS,
and how to transform these parameters to and from other commonly used
triclinic representations.
For style {mesh/hex}, a hexahedral mesh can be read from an ASCII VTK file.
You can apply offset, scaling and rotation to the imported mesh via dedicated keywords.
If applying more then one of these operations, the offset is applied first
and then the geometry is scaled. Then the geometry is rotated around the
x-axis first, then around the y-axis, then around the z-axis.
For style {mesh/tet}, a tetrahedral mesh can be read from an ASCII VTK file.
You can apply offset, scaling and rotation to the imported mesh via dedicated keywords.
If applying more then one of these operations, the offset is applied first

View File

@ -49,12 +49,13 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
ramp(x,y), stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
group functions = count(group), mass(group), charge(group),
xcm(group,dim), vcm(group,dim), fcm(group,dim),
bound(group,xmin), idbound(group,xmin), gyration(group), ke(group),
angmom(group,dim), torque(group,dim),
bound(group,dim), boundid(group,dim), idbound(group,dim), gyration(group),
ke(group), angmom(group,dim), torque(group,dim),
inertia(group,dimdim), omega(group,dim)
region functions = count(group,region), mass(group,region), charge(group,region),
xcm(group,dim,region), vcm(group,dim,region), fcm(group,dim,region),
bound(group,xmin,region), idbound(group,xmin,region), gyration(group,region), ke(group,reigon),
bound(group,dim,region), boundid(group,dim,region), idbound(group,dim,region),
gyration(group,region), ke(group,reigon),
angmom(group,dim,region), torque(group,dim,region),
inertia(group,dimdim,region), omega(group,dim,region)
special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y), next(x)
@ -348,12 +349,12 @@ Math functions: sqrt(x), exp(x), ln(x), log(x), abs(x), \
stagger(x,y), logfreq(x,y,z), stride(x,y,z), vdisplace(x,y), \
swiggle(x,y,z), cwiggle(x,y,z)
Group functions: count(ID), mass(ID), charge(ID), xcm(ID,dim), \
vcm(ID,dim), fcm(ID,dim), bound(ID,dir), idbound(ID,dir),\
gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), \
vcm(ID,dim), fcm(ID,dim), bound(ID,dir), boundid(ID,dir),\
idbound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), \
inertia(ID,dimdim), omega(ID,dim)
Region functions: count(ID,IDR), mass(ID,IDR), charge(ID,IDR), \
xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), \
bound(ID,dir,IDR), idbound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), \
bound(ID,dir,IDR), boundid(ID,dir,IDR), idbound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), \
angmom(ID,dim,IDR), torque(ID,dim,IDR), \
inertia(ID,dimdim,IDR), omega(ID,dim,IDR)
Special functions: sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), \
@ -575,8 +576,9 @@ group functions mass() and charge() are the total mass and charge of
the group. Xcm() and vcm() return components of the position and
velocity of the center of mass of the group. Fcm() returns a
component of the total force on the group of atoms. Bound() returns
the min/max of a particular coordinate for all atoms in the group. Similarly,
IDbound() returns the min/max atom ID in the group.
the min/max of a particular coordinate for all atoms in the group. BoundID()
gives the ID of the particle with the smallest/largest x/y/z coordinate.
Similarly, IDbound() returns the min/max atom ID in the group.
Gyration() computes the radius-of-gyration of the group of atoms. See
the "compute gyration"_compute_gyration.html command for a definition
of the formula. Angmom() returns components of the angular momentum

View File

@ -0,0 +1,158 @@
# 3D bin - fill (coupled coarse grain - original size)
# test for establishing boundary conditions via stress controller
variable cg equal 2
variable cg3 equal ${cg}*${cg}*${cg}
variable rp equal 0.001 # original particle radius
variable dp_mm equal round(2000*${rp})
variable rp_cg equal ${cg}*${rp} # coarse grain radius
variable np_in_reg equal round(360/${cg3}) # number of particles
atom_style sphere
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units si
region reg block -0.021 0.071 -0.021 0.001 -0.021 0.125 units box
create_box 3 reg
neighbor 0.00025 bin
neigh_modify delay 0 exclude type 2 3 # no interaction between particles of type 2 and 3
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8 1.e8 1.e8 # wall coarse resolved
fix m2 all property/global poissonsRatio peratomtype 0.35 0.35 0.35
fix m3 all property/global coefficientRestitution peratomtypepair 3 0.6 0.6 0.6 &
0.6 0.6 0.6 &
0.6 0.6 0.6
fix m4 all property/global coefficientFriction peratomtypepair 3 0.5 0.5 0.5 &
0.5 0.3 0.3 &
0.5 0.3 0.3
fix m5 all property/global coefficientRollingFriction peratomtypepair 3 0.02 0.02 0.02 &
0.02 0.01 0.01 &
0.02 0.01 0.01
# define groups by type
group coarse type 2
group resolved type 3
# pair style
pair_style gran model hertz tangential history rolling_friction epsd2
pair_coeff * *
timestep 0.000001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
# x-walls for coarse part
fix xwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane -0.02
fix xwall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.02
# x-walls for resolved part
fix xwall3 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.03
fix xwall4 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.07
# front & back walls
fix ywall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane -0.02
fix ywall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane 0.00
# bottom wall
fix zwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 zplane -0.02
# generate surface with face ids and the insertion volume, can be removed if files already exist
extrude_surface meshes/surface2x1.vtk file meshes/testsurface2x1.vtk meshes/testextrusion2x1.vtk extrude_length 0.003 min_rad 0.001
extrude_surface meshes/surface2x1_bottom.vtk file meshes/testsurface2x1_bottom.vtk meshes/testextrusion2x1_bottom.vtk extrude_length 0.0035 min_rad 0.0005
# insertion of coarse grain particles
fix pts1 all particletemplate/sphere 1 atom_type 2 density constant 2500 radius constant ${rp_cg}
fix pdd1 all particledistribution/discrete 33335 1 pts1 1.0
region ins_reg block -0.02 0.02 -0.02 0.0 0.10 0.12 units box
fix ins_cg2 coarse insert/pack seed 674562 distributiontemplate pdd1 verbose no &
maxattempt 300 insert_every 500 overlapcheck yes all_in yes vel constant 0. 0. 0. &
region ins_reg particles_in_region ${np_in_reg}
# load surface file with face ids (cell_data option)
fix surface all mesh/surface file meshes/testsurface2x1.vtk type 1 cell_data yes
fix surfaceb all mesh/surface file meshes/testsurface2x1_bottom.vtk type 1 cell_data yes
# measure mass flow through face,
# also specify that all calculations should consider the coarse grain factor (cg 2)
fix massflowcg2 all massflow/mesh/face mesh surface count once cg 2 inside_out #file post/testmassflowprop.txt
fix massflowcg2b all massflow/mesh/face mesh surfaceb count once cg 2 inside_out #file post/testmassflowpropb.txt
# load volume file with face ids for insertion
region hexregion mesh/hex file meshes/testextrusion2x1.vtk scale 1. move 0.05 0. 0. rotate 0. 0. 0. cell_data yes units box
region hexregionb mesh/hex file meshes/testextrusion2x1_bottom.vtk scale 1. move 0.05 0. 0. rotate 0. 0. 0. cell_data yes units box
# insert particles based on the massflow measured,
# also specify that all calculations should consider the coarse grain factor (cg 1)
# and an atom type different from the measured particle shall be used
fix ins_cg1 resolved insert/pack/face seed 7238 random_distribute exact maxattempt 1000 insert_every 500 &
overlapcheck yes all_in yes type_offset 1 region hexregion ntry_mc 10000 massflow_face massflowcg2 cg 1
fix ins_cg1b resolved insert/pack/face seed 72332 random_distribute exact maxattempt 1000 insert_every 500 &
overlapcheck yes all_in yes type_offset 1 region hexregionb ntry_mc 10000 massflow_face massflowcg2b cg 1
# average region of coarse grain simulation
region ave_reg mesh/hex file meshes/grid2x1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
# average region of resolved part
region ave_reg1 mesh/hex file meshes/grid2x1.vtk scale 1. move 0.05 0. 0. rotate 0. 0. 0. cell_data yes units box
variable nevery equal 15
variable one_over_nevery equal 1.0/${nevery}
variable one_over_neverydt equal 1.0/(${nevery}*dt)
# stress computation
fix stress_cg2 coarse ave/euler/region nevery ${nevery} region ave_reg
fix stress_cg1 resolved ave/euler/region nevery ${nevery} region ave_reg1
# control stress in resolved transition region
fix stressctrl resolved forcecontrol/region ctrlPV stress actual_val stress_cg1 target_val stress_cg2 &
kp ${one_over_nevery} ki ${one_over_neverydt} kd 0.0 velocity_limit on cg 2
# deactivate controller in the central part of the resolved region
fix_modify stressctrl deactivate 2 5
# integrator
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
dump dmp0 all custom/vtk 10000 post/bin_fill_d${dp_mm}mm_*.vtk id type radius mass x y z vx vy vz fx fy fz omegax omegay omegaz
#dump dmp1 all euler/vtk 2000 post/eulerresolved*.vtk ave_euler stress_cg1
#dump dmp2 all euler/vtk 2000 post/eulercoarse*.vtk ave_euler stress_cg2
run 140000 upto
# remove any particles above resolved region
region removetop block 0.029 0.071 -0.021 0.001 0.081 0.125 units box
fix removetop_cg1 resolved remove nevery 500 massrate 1000 style delete seed 5143 region removetop atomtype 3
# remove any particles below resolved region
region removebottom block 0.029 0.071 -0.021 0.001 -0.021 -0.001 units box
fix removebottom_cg1 resolved remove nevery 500 massrate 1000 style delete seed 5143 region removebottom atomtype 3
run 1000000 upto
unfix ins_cg1b
unfix massflowcg2b
unfix removebottom_cg1
run 1500000 upto
unfix ins_cg1
unfix massflowcg2
unfix removetop_cg1
run 1760000 upto
unfix ins_cg2
run 2000000 upto

View File

@ -0,0 +1,81 @@
# 3D bin - fill (coarse grain)
variable rp equal 0.002
variable dp_mm equal round(2000*${rp})
variable np_in_reg equal 48
atom_style sphere
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units si
region reg block -0.021 0.021 -0.021 0.001 -0.021 0.125 units box
create_box 2 reg
neighbor 0.00025 bin
neigh_modify delay 0
# material properties required for granular pair style
fix m1 all property/global youngsModulus peratomtype 1.e8 1.e8 # wall particles
fix m2 all property/global poissonsRatio peratomtype 0.35 0.35
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.6 0.6 &
0.6 0.6
fix m4 all property/global coefficientFriction peratomtypepair 2 0.5 0.5 &
0.5 0.3
fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0.02 0.02 &
0.02 0.01
# pair style
pair_style gran model hertz tangential history rolling_friction epsd2
pair_coeff * *
timestep 0.000001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix xwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane -0.02
fix xwall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.02
fix ywall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane -0.02
fix ywall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane 0.0
fix zwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 zplane -0.02
region ins_reg block -0.02 0.02 -0.02 0.0 0.10 0.12 units box
fix pts1 all particletemplate/sphere 1 atom_type 2 density constant 2500 radius constant ${rp}
fix pdd1 all particledistribution/discrete 33335 1 pts1 1.0
fix ins_cg all insert/pack seed 5331 distributiontemplate pdd1 verbose no &
maxattempt 300 insert_every 500 overlapcheck yes all_in yes vel constant 0. 0. 0. &
region ins_reg particles_in_region ${np_in_reg}
# stress computation
#region ave_reg mesh/hex file meshes/grid2x1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
#fix stress_cg all ave/euler/region nevery 10 region ave_reg
# integrator
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
dump dmp0 all custom/vtk 10000 post/bin_fill_d${dp_mm}mm_*.vtk id type radius mass x y z vx vy vz fx fy fz omegax omegay omegaz
#dump dmp1 all euler/vtk 10000 post/euler_cg*.vtk ave_euler stress_cg
run 1760000 upto
unfix ins_cg
run 2000000 upto

View File

@ -0,0 +1,81 @@
# 3D bin - fill (reference size)
variable rp equal 0.001
variable dp_mm equal round(2000*${rp})
variable np_in_reg equal 360
atom_style sphere
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units si
region reg block -0.021 0.021 -0.021 0.001 -0.021 0.125 units box
create_box 2 reg
neighbor 0.00025 bin
neigh_modify delay 0
# material properties required for granular pair style
fix m1 all property/global youngsModulus peratomtype 1.e8 1.e8 # wall particles
fix m2 all property/global poissonsRatio peratomtype 0.35 0.35
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.6 0.6 &
0.6 0.6
fix m4 all property/global coefficientFriction peratomtypepair 2 0.5 0.5 &
0.5 0.3
fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0.02 0.02 &
0.02 0.01
# pair style
pair_style gran model hertz tangential history rolling_friction epsd2
pair_coeff * *
timestep 0.000001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix xwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane -0.02
fix xwall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.02
fix ywall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane -0.02
fix ywall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane 0.0
fix zwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 zplane -0.02
region ins_reg block -0.02 0.02 -0.02 0.0 0.10 0.12 units box
fix pts1 all particletemplate/sphere 1 atom_type 2 density constant 2500 radius constant ${rp}
fix pdd1 all particledistribution/discrete 33335 1 pts1 1.0
fix ins_ref all insert/pack seed 5331 distributiontemplate pdd1 verbose no &
maxattempt 300 insert_every 500 overlapcheck yes all_in yes vel constant 0. 0. 0. &
region ins_reg particles_in_region ${np_in_reg}
# stress computation
#region ave_reg mesh/hex file meshes/grid2x1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
#fix stress_ref all ave/euler/region nevery 10 region ave_reg
# integrator
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
dump dmp0 all custom/vtk 10000 post/bin_fill_d${dp_mm}mm_*.vtk id type radius mass x y z vx vy vz fx fy fz omegax omegay omegaz
#dump dmp1 all euler/vtk 10000 post/euler_ref*.vtk ave_euler stress_ref
run 1760000 upto
unfix ins_ref
run 2000000 upto

View File

@ -0,0 +1,151 @@
# 3D bin - fill (coupled coarse grain - original size)
# test for establishing boundary conditions via stress controller
# run with 2 partitions, e.g.: mpirun -np 3 lmp_ubuntuVTK -in in.bin_fill_universe -partition 1 2
variable rp equal 0.001 # original particle radius
variable dp_mm equal round(2000*${rp})
variable np_in_reg equal 45 # number of cg particles
variable world_name world root child
variable world_cg world 2 1
variable world_send_recv world 2 1
variable world_couplefreq world 500 500
variable world_minx world -0.021 -0.021
variable world_maxx world 0.021 0.021
variable world_miny world -0.021 -0.021
variable world_maxy world 0.001 0.001
variable world_minz world -0.021 -0.001
variable world_maxz world 0.125 0.081
coarsegraining ${world_cg} model_check warn
atom_style sphere
atom_modify map array sort 0 0
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units si
region reg block ${world_minx} ${world_maxx} ${world_miny} ${world_maxy} ${world_minz} ${world_maxz} units box
create_box 3 reg
neighbor 0.00025 bin
neigh_modify delay 0
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8 1.e8 1.e8 # wall coarse resolved
fix m2 all property/global poissonsRatio peratomtype 0.35 0.35 0.35
fix m3 all property/global coefficientRestitution peratomtypepair 3 0.6 0.6 0.6 &
0.6 0.6 0.6 &
0.6 0.6 0.6
fix m4 all property/global coefficientFriction peratomtypepair 3 0.5 0.5 0.5 &
0.5 0.3 0.3 &
0.5 0.3 0.3
fix m5 all property/global coefficientRollingFriction peratomtypepair 3 0.02 0.02 0.02 &
0.02 0.01 0.01 &
0.02 0.01 0.01
# pair style
pair_style gran model hertz tangential history rolling_friction epsd2
pair_coeff * *
timestep 0.000001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
# x-walls for coarse part
fix xwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane -0.02
fix xwall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 xplane 0.02
# front & back walls
fix ywall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane -0.02
fix ywall2 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 yplane 0.00
# bottom wall
fix zwall1 all wall/gran model hertz tangential history rolling_friction epsd2 primitive type 1 zplane -0.02
# generate surface with face ids and the insertion volume, can be removed if files already exist
extrude_surface meshes/surface2x1.vtk file meshes/testsurface2x1.vtk meshes/testextrusion2x1.vtk extrude_length 0.003 min_rad 0.001
extrude_surface meshes/surface2x1_bottom.vtk file meshes/testsurface2x1_bottom.vtk meshes/testextrusion2x1_bottom.vtk extrude_length 0.0035 min_rad 0.0005
# insertion of coarse grain particles
fix pts1 all particletemplate/sphere 1 atom_type 2 density constant 2500 radius constant ${rp}
fix pdd1 all particledistribution/discrete 33335 1 pts1 1.0
partition yes 1 region ins_reg block -0.02 0.02 -0.02 0.0 0.10 0.12 units box
partition yes 1 fix ins_cg2 all insert/pack seed 674562 distributiontemplate pdd1 verbose no &
maxattempt 300 insert_every 500 overlapcheck yes all_in yes vel constant 0. 0. 0. &
region ins_reg particles_in_region ${np_in_reg}
# load surface file with face ids (cell_data option)
partition yes 1 fix surface all mesh/surface file meshes/testsurface2x1.vtk type 1 cell_data yes
partition yes 1 fix surfaceb all mesh/surface file meshes/testsurface2x1_bottom.vtk type 1 cell_data yes
# measure mass flow through face,
# also specify that all calculations should consider the coarse grain factor (cg 2)
partition yes 1 fix massflowcg2 all massflow/mesh/face/universe mesh surface count once inside_out &
send_to_partition ${world_send_recv} couple_every ${world_couplefreq}#file post/testmassflowprop.txt
partition yes 1 fix massflowcg2b all massflow/mesh/face/universe mesh surfaceb count once inside_out &
send_to_partition ${world_send_recv} couple_every ${world_couplefreq}#file post/testmassflowpropb.txt
# load volume file with face ids for insertion
partition yes 2 region hexregion mesh/hex file meshes/testextrusion2x1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
partition yes 2 region hexregionb mesh/hex file meshes/testextrusion2x1_bottom.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
# insert particles based on the massflow measured,
# also specify that all calculations should consider the coarse grain factor (cg 1)
# and an atom type different from the measured particle shall be used
partition yes 2 fix ins_cg1 all insert/pack/face/universe seed 7238 random_distribute exact maxattempt 1000 insert_every ${world_couplefreq} &
overlapcheck yes all_in yes region hexregion ntry_mc 10000 massflow_face massflowcg2 receive_from_partition ${world_send_recv}
partition yes 2 fix ins_cg1b all insert/pack/face/universe seed 72332 random_distribute exact maxattempt 1000 insert_every ${world_couplefreq} &
overlapcheck yes all_in yes region hexregionb ntry_mc 10000 massflow_face massflowcg2b receive_from_partition ${world_send_recv}
# average region of simulation
region ave_reg mesh/hex file meshes/grid2x1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
variable nevery equal 15
variable one_over_nevery equal 1.0/${nevery}
variable one_over_neverydt equal 1.0/(${nevery}*dt)
# stress computation
partition yes 1 fix stress_cg2 all ave/euler/region/universe nevery ${nevery} region ave_reg send_to_partition ${world_send_recv} sync yes
partition yes 2 fix stress_cg1 all ave/euler/region nevery ${nevery} region ave_reg
# control stress in resolved transition region
partition yes 2 fix stressctrl all forcecontrol/region/universe ctrlPV stress actual_val stress_cg1 target_val stress_cg2 &
kp ${one_over_nevery} ki ${one_over_neverydt} kd 0.0 velocity_limit on receive_from_partition ${world_send_recv} couple_every ${nevery}
# deactivate controller in the central part of the resolved region
partition yes 2 fix_modify stressctrl deactivate 2 5
# integrator
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
dump dmp0 all custom/vtk 10000 post/bin_fill_d${dp_mm}mm_${world_name}_*.vtk id type radius mass x y z vx vy vz fx fy fz omegax omegay omegaz
#partition yes 2 dump dmp1 all euler/vtk 2000 post/euler_${world_name}*.vtk ave_euler stress_cg1
#partition yes 1 dump dmp2 all euler/vtk 2000 post/euler_${world_name}*.vtk ave_euler stress_cg2
run 1000000 upto
partition yes 2 unfix ins_cg1b
partition yes 1 unfix massflowcg2b
run 1500000 upto
partition yes 2 unfix ins_cg1
partition yes 1 unfix massflowcg2
run 1760000 upto
partition yes 1 unfix ins_cg2
run 2000000 upto

View File

@ -0,0 +1,77 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 30 float
-0.02 -0.02 0.0
0.0 -0.02 0.0
0.02 -0.02 0.0
-0.02 0.0 0.0
0.0 0.0 0.0
0.02 0.0 0.0
-0.02 -0.02 0.02
0.0 -0.02 0.02
0.02 -0.02 0.02
-0.02 0.0 0.02
0.0 0.0 0.02
0.02 0.0 0.02
-0.02 -0.02 0.04
0.0 -0.02 0.04
0.02 -0.02 0.04
-0.02 0.0 0.04
0.0 0.0 0.04
0.02 0.0 0.04
-0.02 -0.02 0.06
0.0 -0.02 0.06
0.02 -0.02 0.06
-0.02 0.0 0.06
0.0 0.0 0.06
0.02 0.0 0.06
-0.02 -0.02 0.08
0.0 -0.02 0.08
0.02 -0.02 0.08
-0.02 0.0 0.08
0.0 0.0 0.08
0.02 0.0 0.08
CELLS 8 72
8 0 1 4 3 6 7 10 9
8 1 2 5 4 7 8 11 10
8 6 7 10 9 12 13 16 15
8 7 8 11 10 13 14 17 16
8 12 13 16 15 18 19 22 21
8 13 14 17 16 19 20 23 22
8 18 19 22 21 24 25 28 27
8 19 20 23 22 25 26 29 28
CELL_TYPES 8
12
12
12
12
12
12
12
12
CELL_DATA 8
SCALARS cell_id int 1
LOOKUP_TABLE default
0
1
2
3
4
5
6
7
VECTORS stress_ctrl_dir double
0. 0. 1.
0. 0. 1.
0. 0. 0.
0. 0. 0.
0. 0. 0.
0. 0. 0.
0. 0. -1.
0. 0. -1.

View File

@ -0,0 +1,26 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 6 float
-0.02 -0.02 0.08
0.00 -0.02 0.08
0.02 -0.02 0.08
-0.02 0.0 0.08
0.00 0.0 0.08
0.02 0.0 0.08
CELLS 2 10
4 0 1 4 3
4 1 2 5 4
CELL_TYPES 2
9
9
CELL_DATA 2
SCALARS face_id int 1
LOOKUP_TABLE default
6
7

View File

@ -0,0 +1,26 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 6 float
-0.02 -0.02 0.0
0.00 -0.02 0.0
0.02 -0.02 0.0
-0.02 0.0 0.0
0.00 0.0 0.0
0.02 0.0 0.0
CELLS 2 10
4 0 3 4 1
4 1 4 5 2
CELL_TYPES 2
9
9
CELL_DATA 2
SCALARS face_id int 1
LOOKUP_TABLE default
0
1

View File

@ -0,0 +1,121 @@
# hourglass (coupled coarse grain - original size)
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units cgs
region reg block -3.1 3.1 -1.6 1.6 -3.61 3.61 units box
create_box 1 reg
neighbor 0.008 bin
neigh_modify delay 0
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1
region reg_coarse block -3.1 0.0 -1.6 1.6 -3.61 3.61 units box
region reg_resolved block 0.0 3.1 -1.6 1.6 -3.61 3.61 units box
group coarse region reg_coarse
group resolved region reg_resolved
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 981 vector 0.0 0.0 -1.0
variable cg equal 2.0
variable d equal 0.07
variable d_cg equal ${cg}*${d}
variable r equal 0.5*${d}
variable r_cg equal ${cg}*${r}
extract_surface meshes/hourglass_grid.vtk file meshes/hourglass_surface.vtk meshes/hourglass_extrusion.vtk extrude_length 0.19 min_rad ${r}
fix inface all mesh/surface file meshes/hourglass_insert_top.stl type 1 scale 0.9 move -1.55 0.0 0.0
fix hourglass_surf_cg all mesh/surface file meshes/hourglass_surface.vtk type 1 cell_data yes scale 1.0 move -1.55 0.0 0.0
fix hourglass all mesh/surface file meshes/hourglass_wall.stl type 1 scale 1.0 move 1.55 0.0 0.0
fix hourglass_cg all mesh/surface file meshes/hourglass_wall.stl type 1 scale 1.0 move -1.55 0.0 0.0
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 2 meshes hourglass hourglass_cg
fix pts1_cg all particletemplate/sphere 3451 atom_type 1 density constant 2.500 radius constant ${r_cg}
fix pdd1_cg all particledistribution/discrete 5331 1 pts1_cg 1.0
# particle insertion
fix ins coarse insert/stream seed 5331 distributiontemplate pdd1_cg &
nparticles 1500 massrate 100 overlapcheck yes all_in no vel constant 0.0 0.0 -5.0 &
insertion_face inface extrude_length 0.5
region remove_reg_top block 0.005 3.1 -1.6 1.6 2.01 3.6 units box
region remove_reg_bottom block 0.005 3.1 -1.6 1.6 -3.6 -2.01 units box
# measure mass flow through face,
# also specify that all calculations should consider the coarse grain factor (cg 2)
fix massflowcg2 all massflow/mesh/face mesh hourglass_surf_cg count once cg 2 inside_out #file massflow.txt
# load volume file with face ids for insertion
region hexregion mesh/hex file meshes/hourglass_extrusion.vtk scale 1. move 1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
fix ins_cg1 resolved insert/pack/face seed 5330 random_distribute exact maxattempt 150 &
insert_every 500 overlapcheck yes all_in yes type_offset 0 &
region hexregion ntry_mc 10000 massflow_face massflowcg2 cg 1
# average region of coarse grain simulation
region ave_reg mesh/hex file meshes/hourglass_grid.vtk scale 1. move -1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
# average region of resolved part
region ave_reg1 mesh/hex file meshes/hourglass_grid.vtk scale 1. move 1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
variable nevery equal 10
variable one_over_nevery equal 1.0/${nevery}
# stress computation
fix velocity_cg2 coarse ave/euler/region nevery ${nevery} region ave_reg
fix velocity_cg1 resolved ave/euler/region nevery ${nevery} region ave_reg1
# control velocity in coarse region
fix velocityctrl coarse forcecontrol/region ctrlPV velocity actual_val velocity_cg2 target_val velocity_cg1 kp ${one_over_nevery} ki 0.0 kd 0.0 #cg 1
# activate controller only where required
fix_modify velocityctrl deactivate 0 44
fix_modify velocityctrl activate 13
fix_modify velocityctrl activate 22
fix_modify velocityctrl activate 27 35
fix_modify velocityctrl massflow_correction_on 22
# apply nve integration
fix integr all nve/sphere
# output settings
thermo_style custom step atoms ke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 2000 post/hourglass_2lvl*.vtk id type radius x y z vx vy vz fx fy fz
fix removebottom resolved remove nevery 250 massrate 2000 style delete seed 5143 region remove_reg_bottom atomtype 1
run 40000 upto
unfix ins_cg1
unfix massflowcg2
run 400000 upto

View File

@ -0,0 +1,124 @@
# hourglass (coupled coarse grain - original size)
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
processors * 1 1
units cgs
region reg block -3.1 3.1 -1.6 1.6 -3.61 3.61 units box
create_box 1 reg
neighbor 0.008 bin
neigh_modify delay 0
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1
region reg_coarse block -3.1 0.0 -1.6 1.6 -3.61 3.61 units box
region reg_resolved block 0.0 3.1 -1.6 1.6 -3.61 3.61 units box
group coarse region reg_coarse
group resolved region reg_resolved
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 981 vector 0.0 0.0 -1.0
variable cg equal 2.0
variable d equal 0.07
variable d_cg equal ${cg}*${d}
variable r equal 0.5*${d}
variable r_cg equal ${cg}*${r}
extract_surface meshes/hourglass_grid.vtk file meshes/hourglass_surface.vtk meshes/hourglass_extrusion.vtk extrude_length 0.19 min_rad ${r}
fix inface all mesh/surface file meshes/hourglass_insert_top.stl type 1 scale 0.9 move -1.55 0.0 0.0
fix hourglass_surf_cg all mesh/surface file meshes/hourglass_surface.vtk type 1 cell_data yes scale 1.0 move -1.55 0.0 0.0
fix hourglass all mesh/surface file meshes/hourglass_wall.stl type 1 scale 1.0 move 1.55 0.0 0.0
fix hourglass_cg all mesh/surface file meshes/hourglass_wall.stl type 1 scale 1.0 move -1.55 0.0 0.0
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 2 meshes hourglass hourglass_cg
fix pts1_cg all particletemplate/sphere 3451 atom_type 1 density constant 2.500 radius constant ${r_cg}
fix pdd1_cg all particledistribution/discrete 5331 1 pts1_cg 1.0
# particle insertion
fix ins coarse insert/stream seed 5331 distributiontemplate pdd1_cg &
nparticles 1500 massrate 100 overlapcheck yes all_in no vel constant 0.0 0.0 -5.0 &
insertion_face inface extrude_length 0.5
region remove_reg_top block 0.005 3.1 -1.6 1.6 2.01 3.6 units box
region remove_reg_bottom block 0.005 3.1 -1.6 1.6 -3.6 -2.01 units box
# measure mass flow through face,
# also specify that all calculations should consider the coarse grain factor (cg 2)
fix massflowcg2 all massflow/mesh/face mesh hourglass_surf_cg count once cg 2 inside_out #file massflow.txt
# load volume file with face ids for insertion
region hexregion mesh/hex file meshes/hourglass_extrusion.vtk scale 1. move 1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
fix ins_cg1 resolved insert/pack/face seed 5330 random_distribute exact maxattempt 150 &
insert_every 500 overlapcheck yes all_in yes type_offset 0 &
region hexregion ntry_mc 10000 massflow_face massflowcg2 cg 1
# average region of coarse grain simulation
region ave_reg mesh/hex file meshes/hourglass_grid.vtk scale 1. move -1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
# average region of resolved part
region ave_reg1 mesh/hex file meshes/hourglass_grid.vtk scale 1. move 1.55 0. 0. rotate 0. 0. 0. cell_data yes units box
variable nevery equal 10
variable one_over_nevery equal 1.0/${nevery}
# stress computation
fix velocity_cg2 coarse ave/euler/region nevery ${nevery} region ave_reg
fix velocity_cg1 resolved ave/euler/region nevery ${nevery} region ave_reg1
# control velocity in coarse region
fix velocityctrl coarse forcecontrol/region ctrlPV velocity actual_val velocity_cg2 target_val velocity_cg1 kp ${one_over_nevery} ki 0.0 kd 0.0 cg 2
# activate controller only where required
fix_modify velocityctrl deactivate 0 44
fix_modify velocityctrl activate 13
fix_modify velocityctrl activate 22
fix_modify velocityctrl activate 27 35
# correct mass flow by scaling particles inside a spherical region
region scale_reg sphere -1.55 0.0 0.0 0.5 units box
fix scale all scale/diameter 10 region scale_reg scale 1.0 0.7
fix_modify velocityctrl massflow_correction_scale scale 22
# apply nve integration
fix integr all nve/sphere
# output settings
thermo_style custom step atoms ke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 2000 post/hourglass_2lvl_scale*.vtk id type radius x y z vx vy vz fx fy fz
fix removebottom resolved remove nevery 250 massrate 2000 style delete seed 5143 region remove_reg_bottom atomtype 1
run 40000 upto
unfix ins_cg1
unfix massflowcg2
run 400000 upto

View File

@ -0,0 +1,63 @@
# hourglass (coarse grain)
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
processors * * 1
coarsegraining 2.0
units cgs
region reg block -1.6 1.6 -1.6 1.6 -3.61 3.61 units box
create_box 1 reg
neighbor 0.003 bin
neigh_modify delay 0
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 981 vector 0.0 0.0 -1.0
variable d equal 0.07
variable r equal 0.5*${d}
fix hourglasswall all mesh/surface file meshes/hourglass_wall.stl type 1
fix inface all mesh/surface file meshes/hourglass_insert_top.stl type 1 scale 0.9
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes hourglasswall
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2.500 radius constant ${r}
fix pdd1 all particledistribution/discrete 5331 1 pts1 1.0
# particle insertion
fix ins all insert/stream seed 5330 distributiontemplate pdd1 &
nparticles 1500 massrate 100 overlapcheck yes all_in no vel constant 0.0 0.0 -5.0 &
insertion_face inface extrude_length 0.5
# apply nve integration
fix integr all nve/sphere
# output settings
thermo_style custom step atoms ke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 2000 post/hourglass_cg*.vtk id type radius x y z vx vy vz fx fy fz
run 250000 upto

View File

@ -0,0 +1,62 @@
# hourglass (reference size)
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
processors * * 1
units cgs
region reg block -1.6 1.6 -1.6 1.6 -3.61 3.61 units box
create_box 1 reg
neighbor 0.003 bin
neigh_modify delay 0
# material properties required for pair style
fix m1 all property/global youngsModulus peratomtype 1.e8
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 981 vector 0.0 0.0 -1.0
variable d equal 0.07
variable r equal 0.5*${d}
fix hourglasswall all mesh/surface file meshes/hourglass_wall.stl type 1
fix inface all mesh/surface file meshes/hourglass_insert_top.stl type 1 scale 0.9
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes hourglasswall
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2.500 radius constant ${r}
fix pdd1 all particledistribution/discrete 5331 1 pts1 1.0
# particle insertion
fix ins all insert/stream seed 5330 distributiontemplate pdd1 &
nparticles 12000 massrate 100 overlapcheck yes all_in no vel constant 0.0 0.0 -5.0 &
insertion_face inface extrude_length 0.5
# apply nve integration
fix integr all nve/sphere
# output settings
thermo_style custom step atoms ke
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 2000 post/hourglass*.vtk id type radius x y z vx vy vz fx fy fz
run 400000 upto

View File

@ -0,0 +1,245 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 96 float
-1.0 -1.0 -2.0
-0.5 -1.5 -2.0
0.5 -1.5 -2.0
1.0 -1.0 -2.0
-1.5 -0.5 -2.0
-0.5 -0.5 -2.0
0.5 -0.5 -2.0
1.5 -0.5 -2.0
-1.5 0.5 -2.0
-0.5 0.5 -2.0
0.5 0.5 -2.0
1.5 0.5 -2.0
-1.0 1.0 -2.0
-0.5 1.5 -2.0
0.5 1.5 -2.0
1.0 1.0 -2.0
-1.0 -1.0 -1.2
-0.5 -1.5 -1.2
0.5 -1.5 -1.2
1.0 -1.0 -1.2
-1.5 -0.5 -1.2
-0.5 -0.5 -1.2
0.5 -0.5 -1.2
1.5 -0.5 -1.2
-1.5 0.5 -1.2
-0.5 0.5 -1.2
0.5 0.5 -1.2
1.5 0.5 -1.2
-1.0 1.0 -1.2
-0.5 1.5 -1.2
0.5 1.5 -1.2
1.0 1.0 -1.2
-1.0 -1.0 -0.4
-0.5 -1.5 -0.4
0.5 -1.5 -0.4
1.0 -1.0 -0.4
-1.5 -0.5 -0.4
-0.5 -0.5 -0.4
0.5 -0.5 -0.4
1.5 -0.5 -0.4
-1.5 0.5 -0.4
-0.5 0.5 -0.4
0.5 0.5 -0.4
1.5 0.5 -0.4
-1.0 1.0 -0.4
-0.5 1.5 -0.4
0.5 1.5 -0.4
1.0 1.0 -0.4
-1.0 -1.0 0.4
-0.5 -1.5 0.4
0.5 -1.5 0.4
1.0 -1.0 0.4
-1.5 -0.5 0.4
-0.5 -0.5 0.4
0.5 -0.5 0.4
1.5 -0.5 0.4
-1.5 0.5 0.4
-0.5 0.5 0.4
0.5 0.5 0.4
1.5 0.5 0.4
-1.0 1.0 0.4
-0.5 1.5 0.4
0.5 1.5 0.4
1.0 1.0 0.4
-1.0 -1.0 1.2
-0.5 -1.5 1.2
0.5 -1.5 1.2
1.0 -1.0 1.2
-1.5 -0.5 1.2
-0.5 -0.5 1.2
0.5 -0.5 1.2
1.5 -0.5 1.2
-1.5 0.5 1.2
-0.5 0.5 1.2
0.5 0.5 1.2
1.5 0.5 1.2
-1.0 1.0 1.2
-0.5 1.5 1.2
0.5 1.5 1.2
1.0 1.0 1.2
-1.0 -1.0 2.0
-0.5 -1.5 2.0
0.5 -1.5 2.0
1.0 -1.0 2.0
-1.5 -0.5 2.0
-0.5 -0.5 2.0
0.5 -0.5 2.0
1.5 -0.5 2.0
-1.5 0.5 2.0
-0.5 0.5 2.0
0.5 0.5 2.0
1.5 0.5 2.0
-1.0 1.0 2.0
-0.5 1.5 2.0
0.5 1.5 2.0
1.0 1.0 2.0
CELLS 45 405
8 0 1 5 4 16 17 21 20
8 1 2 6 5 17 18 22 21
8 2 3 7 6 18 19 23 22
8 4 5 9 8 20 21 25 24
8 5 6 10 9 21 22 26 25
8 6 7 11 10 22 23 27 26
8 8 9 13 12 24 25 29 28
8 9 10 14 13 25 26 30 29
8 10 11 15 14 26 27 31 30
8 16 17 21 20 32 33 37 36
8 17 18 22 21 33 34 38 37
8 18 19 23 22 34 35 39 38
8 20 21 25 24 36 37 41 40
8 21 22 26 25 37 38 42 41
8 22 23 27 26 38 39 43 42
8 24 25 29 28 40 41 45 44
8 25 26 30 29 41 42 46 45
8 26 27 31 30 42 43 47 46
8 32 33 37 36 48 49 53 52
8 33 34 38 37 49 50 54 53
8 34 35 39 38 50 51 55 54
8 36 37 41 40 52 53 57 56
8 37 38 42 41 53 54 58 57
8 38 39 43 42 54 55 59 58
8 40 41 45 44 56 57 61 60
8 41 42 46 45 57 58 62 61
8 42 43 47 46 58 59 63 62
8 48 49 53 52 64 65 69 68
8 49 50 54 53 65 66 70 69
8 50 51 55 54 66 67 71 70
8 52 53 57 56 68 69 73 72
8 53 54 58 57 69 70 74 73
8 54 55 59 58 70 71 75 74
8 56 57 61 60 72 73 77 76
8 57 58 62 61 73 74 78 77
8 58 59 63 62 74 75 79 78
8 64 65 69 68 80 81 85 84
8 65 66 70 69 81 82 86 85
8 66 67 71 70 82 83 87 86
8 68 69 73 72 84 85 89 88
8 69 70 74 73 85 86 90 89
8 70 71 75 74 86 87 91 90
8 72 73 77 76 88 89 93 92
8 73 74 78 77 89 90 94 93
8 74 75 79 78 90 91 95 94
CELL_TYPES 45
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
CELL_DATA 45
SCALARS cell_id int 1
LOOKUP_TABLE default
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

View File

@ -0,0 +1,128 @@
solid ascii
facet normal 0 0 1
outer loop
vertex -0.5 -1.5 3.1
vertex -0.5 -0.5 3.1
vertex -1 -1 3.1
endloop
endfacet
facet normal 0 -0 1
outer loop
vertex -1.5 -0.5 3.1
vertex -1 -1 3.1
vertex -0.5 -0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex -0.5 -1.5 3.1
vertex 0.5 -1.5 3.1
vertex -0.5 -0.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex 0.5 -0.5 3.1
vertex -0.5 -0.5 3.1
vertex 0.5 -1.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex 0.5 -1.5 3.1
vertex 1 -1 3.1
vertex 0.5 -0.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex 1.5 -0.5 3.1
vertex 0.5 -0.5 3.1
vertex 1 -1 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex -1.5 -0.5 3.1
vertex -0.5 -0.5 3.1
vertex -1.5 0.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex -0.5 0.5 3.1
vertex -1.5 0.5 3.1
vertex -0.5 -0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex -0.5 -0.5 3.1
vertex 0.5 -0.5 3.1
vertex -0.5 0.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex 0.5 0.5 3.1
vertex -0.5 0.5 3.1
vertex 0.5 -0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex 0.5 -0.5 3.1
vertex 1.5 -0.5 3.1
vertex 0.5 0.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex 1.5 0.5 3.1
vertex 0.5 0.5 3.1
vertex 1.5 -0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex -1.5 0.5 3.1
vertex -0.5 0.5 3.1
vertex -1 1 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex -0.5 1.5 3.1
vertex -1 1 3.1
vertex -0.5 0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex -0.5 0.5 3.1
vertex 0.5 0.5 3.1
vertex -0.5 1.5 3.1
endloop
endfacet
facet normal -0 0 1
outer loop
vertex 0.5 1.5 3.1
vertex -0.5 1.5 3.1
vertex 0.5 0.5 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex 0.5 1.5 3.1
vertex 0.5 0.5 3.1
vertex 1 1 3.1
endloop
endfacet
facet normal 0 0 1
outer loop
vertex 1.5 0.5 3.1
vertex 1 1 3.1
vertex 0.5 0.5 3.1
endloop
endfacet
endsolid

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,22 @@
{
"runs": [
{
"name" : "hourglass_ref",
"input_script" : "in.hourglass_ref",
"type" : "mpi",
"nprocs" : 4
},
{
"name" : "hourglass_cg",
"input_script": "in.hourglass_cg",
"type": "mpi",
"nprocs": 4
},
{
"name" : "hourglass_2lvl",
"input_script": "in.hourglass_2lvl",
"type": "mpi",
"nprocs": 4
}
]
}

View File

@ -0,0 +1,72 @@
# ave/euler example
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -0.5 0.5 -0.5 0.5 0.0 1.5 units box
create_box 1 reg
neighbor 0.025 bin
neigh_modify delay 0
# material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e7
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.7
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
# pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.000005
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix xwalls1 all wall/gran model hertz tangential history primitive type 1 xplane -0.5
fix xwalls2 all wall/gran model hertz tangential history primitive type 1 xplane 0.5
fix ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane -0.5
fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane 0.5
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0
# region and insertion
variable r equal 0.04
region ins_reg block -0.5 0.5 -0.5 0.5 1.3 1.5 units box
fix pts1 all particletemplate/sphere 33331 atom_type 1 density constant 2000 radius constant ${r}
fix pdd1 all particledistribution/discrete 5331 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.1 &
insert_every 2000 overlapcheck yes all_in no volumefraction_region 0.3 region ins_reg
# stress computation
fix stress all ave/euler nevery 10 cell_size_relative 6 6 6 parallel no
# apply nve integration
fix integr all nve/sphere
# output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
# insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 2000 post/ave_euler*.vtk id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
dump dmpstress all euler/vtk 2000 post/stress_ave_euler*.vtk
run 250000
unfix ins
run 150000

View File

@ -0,0 +1,73 @@
# ave/euler/region example
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -0.5 0.5 -0.5 0.5 0.0 1.5 units box
create_box 1 reg
neighbor 0.025 bin
neigh_modify delay 0
# material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e7
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.7
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
# pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.000005
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix xwalls1 all wall/gran model hertz tangential history primitive type 1 xplane -0.5
fix xwalls2 all wall/gran model hertz tangential history primitive type 1 xplane 0.5
fix ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane -0.5
fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane 0.5
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0
# region and insertion
variable r equal 0.04
region ave_reg mesh/hex file meshes/test.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
region ins_reg block -0.5 0.5 -0.5 0.5 1.3 1.5 units box
fix pts1 all particletemplate/sphere 33331 atom_type 1 density constant 2000 radius constant ${r}
fix pdd1 all particledistribution/discrete 5331 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.1 &
insert_every 2000 overlapcheck yes all_in no volumefraction_region 0.3 region ins_reg
# stress computation
fix stress all ave/euler/region nevery 10 region ave_reg
# apply nve integration
fix integr all nve/sphere
# output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 2000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
# insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 2000 post/ave_euler_region*.vtk id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
dump dmpstress all euler/vtk 2000 post/stress_ave_euler_region*.vtk
run 250000
unfix ins
run 150000

View File

@ -0,0 +1,85 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 36 float
-0.5 -0.5 0.0
0.0 -0.5 0.0
0.5 -0.5 0.0
-0.5 0.0 0.0
0.0 0.0 0.0
0.5 0.0 0.0
-0.5 0.5 0.0
0.0 0.5 0.0
0.5 0.5 0.0
-0.5 -0.5 0.5
0.0 -0.5 0.5
0.5 -0.5 0.5
-0.5 0.0 0.5
0.0 0.0 0.5
0.5 0.0 0.5
-0.5 0.5 0.5
0.0 0.5 0.5
0.5 0.5 0.5
-0.5 -0.5 1.0
0.0 -0.5 1.0
0.5 -0.5 1.0
-0.5 0.0 1.0
0.0 0.0 1.0
0.5 0.0 1.0
-0.5 0.5 1.0
0.0 0.5 1.0
0.5 0.5 1.0
-0.5 -0.5 1.5
0.0 -0.5 1.5
0.5 -0.5 1.5
-0.5 0.0 1.5
0.0 0.0 1.5
0.5 0.0 1.5
-0.5 0.5 1.5
0.0 0.5 1.5
0.5 0.5 1.5
CELLS 12 108
8 0 1 4 3 9 10 13 12
8 1 2 5 4 10 11 14 13
8 3 4 7 6 12 13 16 15
8 4 5 8 7 13 14 17 16
8 9 10 13 12 18 19 22 21
8 10 11 14 13 19 20 23 22
8 12 13 16 15 21 22 25 24
8 13 14 17 16 22 23 26 25
8 18 19 22 21 27 28 31 30
8 19 20 23 22 28 29 32 31
8 21 22 25 24 30 31 34 33
8 22 23 26 25 31 32 35 34
CELL_TYPES 12
12
12
12
12
12
12
12
12
12
12
12
12
CELL_DATA 12
SCALARS cell_id int 1
LOOKUP_TABLE default
0
1
2
3
4
5
6
7
8
9
10
11

View File

View File

@ -0,0 +1,14 @@
{
"runs": [
{
"name" : "ave_euler",
"input_script" : "in.ave_euler",
"type" : "serial"
},
{
"name" : "ave_euler_region",
"input_script": "in.ave_euler_region",
"type": "serial"
}
]
}

View File

@ -0,0 +1,65 @@
#Neighbor list exclusion example
atom_style granular
atom_modify map array
boundary m m m
newton off
communicate single vel yes
units si
region reg block -0.05 0.05 -0.05 0.05 0. 0.15 units box
create_box 2 reg
neighbor 0.001 bin
neigh_modify delay 0 exclude type 1 2
#material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.05 0.05 0.05 0.05
#pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0
fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 0.15
fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.05 0. 0.
#region of insertion
region bc cylinder z 0. 0. 0.045 0.00 0.15 units box
#particle distributions
fix pts1 all particletemplate/sphere 451 atom_type 1 density constant 2500 radius constant 0.002
fix pts2 all particletemplate/sphere 671 atom_type 2 density constant 2500 radius constant 0.004
fix pdd1 all particledistribution/discrete 45341 2 pts1 0.5 pts2 0.5
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. -0.5 &
insert_every once overlapcheck yes all_in yes particles_in_region 5400 region bc
#apply nve integration to all particles
fix integr all nve/sphere
#output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
dump dmp all custom/vtk 800 post/exclude*.vtk id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
unfix ins
run 50000 upto

View File

View File

@ -0,0 +1,9 @@
{
"runs": [
{
"name" : "exclude_neighs",
"input_script" : "in.exclude",
"type" : "serial"
}
]
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,50 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 27 float
-0.5 -0.5 0.0
0.0 -0.5 0.0
0.5 -0.5 0.0
-0.5 0.0 0.0
0.0 0.0 0.0
0.5 0.0 0.0
-0.5 0.5 0.0
0.0 0.5 0.0
0.5 0.5 0.0
-0.5 -0.5 0.5
0.0 -0.5 0.5
0.5 -0.5 0.5
-0.5 0.0 0.5
0.0 0.0 0.5
0.5 0.0 0.5
-0.5 0.5 0.5
0.0 0.5 0.5
0.5 0.5 0.5
-0.5 -0.5 1.0
0.0 -0.5 1.0
0.5 -0.5 1.0
-0.5 0.0 1.0
0.0 0.0 1.0
0.5 0.0 1.0
-0.5 0.5 1.0
0.0 0.5 1.0
0.5 0.5 1.0
CELLS 8 72
8 0 1 4 3 9 10 13 12
8 1 2 5 4 10 11 14 13
8 3 4 7 6 12 13 16 15
8 4 5 8 7 13 14 17 16
8 9 10 13 12 18 19 22 21
8 10 11 14 13 19 20 23 22
8 12 13 16 15 21 22 25 24
8 13 14 17 16 22 23 26 25
CELL_TYPES 8
12
12
12
12
12
12
12
12

View File

@ -0,0 +1,52 @@
# extract surface from hexhedral mesh example
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.05 1.05 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
#material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
#pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
extract_surface data/test.vtk file data/testsurface.vtk data/testextrusion.vtk extrude_length 0.1 min_rad 0.025
#region and insertion
variable d equal 0.05
create_atoms 1 single 0 0 0.5
set group all density 2000 diameter $d
#apply nve integration
fix integr all nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1

View File

@ -0,0 +1,52 @@
# extract surface from hexhedral mesh example
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.05 1.05 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
#material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
#pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
extrude_surface data/shell.vtk file data/testsurface.vtk data/testextrusion.vtk extrude_length 0.1 min_rad 0.025
#region and insertion
variable d equal 0.05
create_atoms 1 single 0 0 0.5
set group all density 2000 diameter $d
#apply nve integration
#fix integr all nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1 post no

View File

@ -0,0 +1,96 @@
# measure mass flow through mesh surface and use it for insertion of particles with a different resolution
variable my_timestep equal 0.00001
variable my_timesteps equal 60000
variable my_dumpfreq equal 1000
variable my_couplefreq equal 400
variable my_minx equal -1
variable my_maxx equal 1
variable my_miny equal -1
variable my_maxy equal 1
variable my_minz equal -0.5
variable my_maxz equal 1.5
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block ${my_minx} ${my_maxx} ${my_miny} ${my_maxy} ${my_minz} ${my_maxz} units box
create_box 2 reg
neighbor 0.005 bin
# 1) mutually exclude types representing different resolution levels from neighbor list build
neigh_modify delay 0 exclude type 1 2
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.05 0.05 0.05 0.05
# 2) define a property that will hold the coarse grain level for each particle
fix m5 all property/atom/lammps i_cg_level
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep ${my_timestep}
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix pts1 all particletemplate/sphere 98531 atom_type 1 density constant 2000 radius constant 0.02
fix pdd1 all particledistribution/discrete 746325 1 pts1 1.0
# region to insert cg 4 particles
region insregioncg4 block -0.65 0.65 -0.65 0.65 1.1 1.3 units box
# cg 4 particle insertion
fix ins_cg4 all insert/pack seed 5330 distributiontemplate pdd1 &
maxattempt 200 insert_every once overlapcheck yes all_in yes set_property i_cg_level 4 &
vel constant 0. 0. 0. region insregioncg4 volumefraction_region 0.1
# 3) generate surface with face ids and the insertion volume, can be removed if files already exist
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.035 min_rad 0.01
# 4) load surface file with face ids (cell_data option)
fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
# 5) measure mass flow through face restricting to particles with property i_cg_level 4,
# also specify that all calculations should consider the coarse grain factor (cg 4)
fix massflowcg4 all massflow/mesh/face mesh surface count once check_property i_cg_level 4 cg 4 inside_out #file post/testmassflowprop.txt
# 6) load volume file with face ids for insertion
region hexregion mesh/hex file meshes/testextrusion.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
# 7) insert particles based on the massflow measured, set property i_cg_level 2,
# also specify that all calculations should consider the coarse grain factor (cg 2)
# and an atom type different from the measured particle shall be used
fix ins_cg2 all insert/pack/face seed 7331 random_distribute exact maxattempt 250 insert_every ${my_couplefreq} &
overlapcheck yes set_property i_cg_level 2 all_in yes type_offset 1 region hexregion ntry_mc 10000 massflow_face massflowcg4 cg 2
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes extra 0
run 1
dump dmp all custom/vtk ${my_dumpfreq} post/insert*.vtk id type radius mass x y z vx vy vz fx fy fz i_cg_level
region removeregion block -0.5 0.5 -0.5 0.5 0.0 1.0 side out units box
fix remove_cg2 all remove nevery 400 massrate 1000 style delete seed 5143 region removeregion atomtype 2
run ${my_timesteps} upto

View File

@ -0,0 +1,169 @@
# measure mass flow through mesh surface and use it for insertion of particles with a different resolution
# additionally measure temperature and relative layer radii and set them upon insertion
variable my_timestep equal 0.00001
variable my_timesteps equal 60000
variable my_dumpfreq equal 1000
variable my_couplefreq equal 400
variable my_minx equal -1
variable my_maxx equal 1
variable my_miny equal -1
variable my_maxy equal 1
variable my_minz equal -0.5
variable my_maxz equal 1.5
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block ${my_minx} ${my_maxx} ${my_miny} ${my_maxy} ${my_minz} ${my_maxz} units box
create_box 2 reg
neighbor 0.005 bin
# 1) mutually exclude types representing different resolution levels from neighbor list build
neigh_modify delay 0 exclude type 1 2
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.05 0.05 0.05 0.05
# 2) define a property that will hold the coarse grain level for each particle
fix m5 all property/atom/lammps i_cg_level
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep ${my_timestep}
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
group cg4group region reg
group cg2group region reg
################################################################################
# heat transfer
fix ftco all property/global thermalConductivity peratomtype 100. 100.
fix ftca all property/global thermalCapacity peratomtype 10. 10.
fix heattransfer all heat/gran initial_temperature 1100.
################################################################################
# chemistry
# input fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix partTemp all property/atom partTemp scalar yes no no 1223.15 # Kelvin
fix partNu all property/atom partNu scalar yes no no 0.00004125
fix partRe all property/atom partRe scalar yes no no 2.7500
fix partP all property/atom partP scalar yes no no 101325
# molar fractions
fix X_CO all property/atom X_CO scalar yes no no 0.4
fix X_CO2 all property/atom X_CO2 scalar yes no no 0.0
fix X_N2 all property/atom X_N2 scalar yes no no 0.6
# molecular diffusion coefficient
fix CO_diffCoeff all property/atom CO_diffCoeff scalar yes no no 0.00023
# output fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix Modified_CO all property/atom Modified_CO scalar yes yes no 0.0
fix Modified_CO2 all property/atom Modified_CO2 scalar yes yes no 0.0
fix reactionHeat all property/atom reactionHeat scalar yes no no 0.0
################################################################################
# activate 3-layer unreacted shrinking core model
fix cfd54 cg4group chem/shrink/core speciesA CO molMassA 0.02801 &
speciesC CO2 molMassC 0.04401 &
scale_reduction_rate 10.0 nevery 10 cg 4 #screen yes
fix cfd52 cg2group chem/shrink/core speciesA CO molMassA 0.02801 &
speciesC CO2 molMassC 0.04401 &
scale_reduction_rate 10.0 nevery 10 cg 2 #screen yes
# chemical properties for unreacted shrinking core model
# w->Fe m->w h->m
fix k0_CO4 all property/global k0_cfd54 vector 17 25 2700
fix Ea_CO4 all property/global Ea_cfd54 vector 69488 73674 113859
fix k0_CO2 all property/global k0_cfd52 vector 17 25 2700
fix Ea_CO2 all property/global Ea_cfd52 vector 69488 73674 113859
# particle porosity/tortuosity/pore diameter
fix porosity4 all property/global porosity_cg4group vector 0.61 0.31 0.16 0.15
fix tortuosity4 all property/global tortuosity_cg4group scalar 3
fix pore_diameter4 all property/global pore_diameter_cg4group scalar 5.5e-7
fix porosity2 all property/global porosity_cg2group vector 0.61 0.31 0.16 0.15
fix tortuosity2 all property/global tortuosity_cg2group scalar 3
fix pore_diameter2 all property/global pore_diameter_cg2group scalar 5.5e-7
# define layer properties
fix LayerRelRadii all property/atom relRadii vector yes no no 1.0 0.998 0.995 0.98
# define fix for layer densities
fix LayerDensities4 all property/global density_cg4group vector 7870. 5740. 5170. 5240.
fix LayerDensities2 all property/global density_cg2group vector 7870. 5740. 5170. 5240.
################################################################################
fix pts1 all particletemplate/sphere 98531 atom_type 1 density constant 2000 radius constant 0.02
fix pdd1 all particledistribution/discrete 746325 1 pts1 1.0
# region to insert cg 4 particles
region insregioncg4 block -0.65 0.65 -0.65 0.65 1.1 1.3 units box
# cg 4 particle insertion
fix ins_cg4 cg4group insert/pack seed 5330 distributiontemplate pdd1 &
maxattempt 200 insert_every once overlapcheck yes all_in yes set_property i_cg_level 4 &
vel constant 0. 0. 0. region insregioncg4 volumefraction_region 0.1
# 3) generate surface with face ids and the insertion volume, can be removed if files already exist
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.035 min_rad 0.01
# 4) load surface file with face ids (cell_data option)
fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
# 5) measure mass flow through face restricting to particles with property i_cg_level 4,
# also specify that all calculations should consider the coarse grain factor (cg 4)
fix massflowcg4 cg4group massflow/mesh/face mesh surface count once &
check_property i_cg_level 4 cg 4 inside_out temperature yes chemistry yes #file post/testmassflowprop.txt
# 6) load volume file with face ids for insertion
region hexregion mesh/hex file meshes/testextrusion.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
# 7) insert particles based on the massflow measured, set property i_cg_level 2,
# also specify that all calculations should consider the coarse grain factor (cg 2)
# and an atom type different from the measured particle shall be used
fix ins_cg2 cg2group insert/pack/face seed 7331 random_distribute exact maxattempt 250 insert_every ${my_couplefreq} &
overlapcheck yes set_property i_cg_level 2 all_in yes type_offset 1 temperature yes chemistry yes &
region hexregion ntry_mc 10000 massflow_face massflowcg4 cg 2
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes extra 0
run 1
variable fracRedTotal atom 1./9*f_fracRed[1]+2./9.*f_fracRed[2]+6./9.*f_fracRed[3]
dump dmp all custom/vtk ${my_dumpfreq} post/insert*.vtk id type radius mass &
x y z vx vy vz fx fy fz f_Temp f_fracRed[1] f_fracRed[2] f_fracRed[3] &
f_LayerRelRadii[1] f_LayerRelRadii[2] f_LayerRelRadii[3] f_LayerRelRadii[4] &
v_fracRedTotal i_cg_level
region highTregion block -0.65 0.65 -0.65 0.65 1.2 1.3 units box
group hot region highTregion
set group hot property/atom Temp 1200
region removeregion block -0.5 0.5 -0.5 0.5 0.0 1.0 side out units box
fix remove_cg2 all remove nevery 400 massrate 1000 style delete seed 5143 region removeregion atomtype 2
run ${my_timesteps} upto

View File

@ -0,0 +1,150 @@
# measure mass flow through mesh surface and use it for insertion of particles with a different resolution
# run with 2 partitions, e.g.: mpirun -np 4 lmp_ubuntuVTK -in in.insert_pack_face_ex_universe -partition 2 2
variable world_name world root child
variable world_partition world 1 2
variable world_send_recv world 2 1
variable world_timestep world 0.00001 0.00001 # 0.00001 0.000005
variable world_timesteps world 60000 60000 # 60000 120000
variable world_dumpfreq world 1000 1000 #500 1000
variable world_couplefreq world 400 400 #400 800
variable world_cg world 4 2
variable world_minx world -1 -0.51
variable world_maxx world 1 0.51
variable world_miny world -1 -0.51
variable world_maxy world 1 0.51
variable world_minz world -0.5 -0.01
variable world_maxz world 1.5 1.01
coarsegraining ${world_cg} model_check warn
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block ${world_minx} ${world_maxx} ${world_miny} ${world_maxy} ${world_minz} ${world_maxz} units box
create_box 2 reg
neighbor 0.005 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.05 0.05 0.05 0.05
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep ${world_timestep}
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
################################################################################
# heat transfer
fix ftco all property/global thermalConductivity peratomtype 100. 100.
fix ftca all property/global thermalCapacity peratomtype 10. 10.
fix heattransfer all heat/gran initial_temperature 1100.
################################################################################
# chemistry
# input fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix partTemp all property/atom partTemp scalar yes no no 1223.15 # Kelvin
fix partNu all property/atom partNu scalar yes no no 0.00004125
fix partRe all property/atom partRe scalar yes no no 2.7500
fix partP all property/atom partP scalar yes no no 101325
# molar fractions
fix X_CO all property/atom X_CO scalar yes no no 0.4
fix X_CO2 all property/atom X_CO2 scalar yes no no 0.0
fix X_N2 all property/atom X_N2 scalar yes no no 0.6
# molecular diffusion coefficient
fix CO_diffCoeff all property/atom CO_diffCoeff scalar yes no no 0.00023
# output fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix Modified_CO all property/atom Modified_CO scalar yes yes no 0.0
fix Modified_CO2 all property/atom Modified_CO2 scalar yes yes no 0.0
fix reactionHeat all property/atom reactionHeat scalar yes no no 0.0
################################################################################
# activate 3-layer unreacted shrinking core model
fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 &
speciesC CO2 molMassC 0.04401 &
scale_reduction_rate 10.0 nevery 10
# chemical properties for unreacted shrinking core model
# w->Fe m->w h->m
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.61 0.31 0.16 0.15
fix tortuosity all property/global tortuosity_all scalar 3
fix pore_diameter all property/global pore_diameter_all scalar 5.5e-7
# define layer properties
fix LayerRelRadii all property/atom relRadii vector yes no no 1.0 0.998 0.995 0.98
# define fix for layer densities
fix LayerDensities all property/global density_all vector 7870. 5740. 5170. 5240.
################################################################################
fix pts1 all particletemplate/sphere 98531 atom_type 1 density constant 2000 radius constant 0.005
fix pdd1 all particledistribution/discrete 746325 1 pts1 1.0
# 1) generate surface with face ids and the insertion volume, can be removed if files already exist
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.035 min_rad 0.01
# cg 4 (partition 1):
# region to insert particles
# particle insertion
# 2) load surface file with face ids (cell_data option)
# 3) measure mass flow through face, specify partition to send data to
partition yes 1 region insregioncg4 block -0.65 0.65 -0.65 0.65 1.1 1.3 units box
partition yes 1 fix ins_cg4 all insert/pack seed 5330 distributiontemplate pdd1 maxattempt 200 insert_every once overlapcheck yes all_in yes vel constant 0. 0. 0. region insregioncg4 volumefraction_region 0.1
partition yes 1 fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
partition yes 1 fix massflowcg4 all massflow/mesh/face/universe mesh surface count once inside_out &
temperature yes chemistry yes send_to_partition ${world_send_recv} couple_every ${world_couplefreq}
# cg 2 (partition 2):
# 4) load volume file with face ids for insertion
# 5) insert particles based on the massflow measured, specify partition to receive data from
partition yes 2 region hexregion mesh/hex file meshes/testextrusion.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
partition yes 2 fix ins_cg2 all insert/pack/face/universe seed 7331 random_distribute exact maxattempt 250 insert_every ${world_couplefreq} &
overlapcheck yes all_in yes region hexregion ntry_mc 10000 massflow_face massflowcg4 &
temperature yes chemistry yes receive_from_partition ${world_send_recv}
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes extra 0
run 1
variable fracRedTotal atom 1./9*f_fracRed[1]+2./9.*f_fracRed[2]+6./9.*f_fracRed[3]
dump dmp all custom/vtk ${world_dumpfreq} post/insert_${world_name}_*.vtk id type radius mass &
x y z vx vy vz fx fy fz f_Temp f_fracRed[1] f_fracRed[2] f_fracRed[3] &
f_LayerRelRadii[1] f_LayerRelRadii[2] f_LayerRelRadii[3] f_LayerRelRadii[4] &
v_fracRedTotal
partition yes 1 region highTregion block -0.65 0.65 -0.65 0.65 1.2 1.3 units box
partition yes 1 group hot region highTregion
partition yes 1 set group hot property/atom Temp 1200
run ${world_timesteps} upto

View File

@ -0,0 +1,90 @@
# measure mass flow through mesh surface and use it for insertion of particles with a different resolution
# run with 2 partitions, e.g.: mpirun -np 4 lmp_ubuntuVTK -in in.insert_pack_face_universe -partition 2 2
variable world_name world root child
variable world_partition world 1 2
variable world_send_recv world 2 1
variable world_timestep world 0.00001 0.00001 # 0.00001 0.000005
variable world_timesteps world 60000 60000 # 60000 120000
variable world_dumpfreq world 1000 1000 #500 1000
variable world_couplefreq world 400 400 #400 800
variable world_cg world 4 2
variable world_minx world -1 -0.51
variable world_maxx world 1 0.51
variable world_miny world -1 -0.51
variable world_maxy world 1 0.51
variable world_minz world -0.5 -0.01
variable world_maxz world 1.5 1.01
coarsegraining ${world_cg}
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block ${world_minx} ${world_maxx} ${world_miny} ${world_maxy} ${world_minz} ${world_maxz} units box
create_box 2 reg
neighbor 0.005 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.9 0.85 0.85 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.05 0.05 0.05 0.05
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep ${world_timestep}
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
fix pts1 all particletemplate/sphere 98531 atom_type 1 density constant 2000 radius constant 0.005
fix pdd1 all particledistribution/discrete 746325 1 pts1 1.0
# 1) generate surface with face ids and the insertion volume, can be removed if files already exist
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.035 min_rad 0.01
# cg 4 (partition 1):
# region to insert particles
# particle insertion
# 2) load surface file with face ids (cell_data option)
# 3) measure mass flow through face, specify partition to send data to
partition yes 1 region insregioncg4 block -0.65 0.65 -0.65 0.65 1.1 1.3 units box
partition yes 1 fix ins_cg4 all insert/pack seed 5330 distributiontemplate pdd1 maxattempt 200 insert_every once overlapcheck yes all_in yes vel constant 0. 0. 0. region insregioncg4 volumefraction_region 0.1
partition yes 1 fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
partition yes 1 fix massflowcg4 all massflow/mesh/face/universe mesh surface count once inside_out send_to_partition ${world_send_recv} couple_every ${world_couplefreq}
# cg 2 (partition 2):
# 4) load volume file with face ids for insertion
# 5) insert particles based on the massflow measured, specify partition to receive data from
partition yes 2 region hexregion mesh/hex file meshes/testextrusion.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. cell_data yes units box
partition yes 2 fix ins_cg2 all insert/pack/face/universe seed 7331 random_distribute exact maxattempt 250 insert_every ${world_couplefreq} &
overlapcheck yes all_in yes region hexregion ntry_mc 10000 massflow_face massflowcg4 receive_from_partition ${world_send_recv}
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes extra 0
run 1
dump dmp all custom/vtk ${world_dumpfreq} post/insert_${world_name}_*.vtk id type radius mass x y z vx vy vz fx fy fz
run ${world_timesteps} upto

View File

@ -0,0 +1,50 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 27 float
-0.5 -0.5 0.0
0.0 -0.5 0.0
0.5 -0.5 0.0
-0.5 0.0 0.0
0.0 0.0 0.0
0.5 0.0 0.0
-0.5 0.5 0.0
0.0 0.5 0.0
0.5 0.5 0.0
-0.5 -0.5 0.5
0.0 -0.5 0.5
0.5 -0.5 0.5
-0.5 0.0 0.5
0.0 0.0 0.5
0.5 0.0 0.5
-0.5 0.5 0.5
0.0 0.5 0.5
0.5 0.5 0.5
-0.5 -0.5 1.0
0.0 -0.5 1.0
0.5 -0.5 1.0
-0.5 0.0 1.0
0.0 0.0 1.0
0.5 0.0 1.0
-0.5 0.5 1.0
0.0 0.5 1.0
0.5 0.5 1.0
CELLS 8 72
8 0 1 4 3 9 10 13 12
8 1 2 5 4 10 11 14 13
8 3 4 7 6 12 13 16 15
8 4 5 8 7 13 14 17 16
8 9 10 13 12 18 19 22 21
8 10 11 14 13 19 20 23 22
8 12 13 16 15 21 22 25 24
8 13 14 17 16 22 23 26 25
CELL_TYPES 8
12
12
12
12
12
12
12
12

View File

@ -0,0 +1,14 @@
{
"runs": [
{
"name" : "insert_pack_face",
"input_script" : "in.insert_pack_face",
"type" : "serial"
},
{
"name" : "insert_pack_face_ex",
"input_script": "in.insert_pack_face_ex",
"type": "serial"
}
]
}

View File

@ -0,0 +1,219 @@
# measure mass flow through mesh surface (broken down to individual faces)
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.2 1.2 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.8
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 1.0 1.0 -1.0
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.1 min_rad 0.025
fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
fix massflowout all massflow/mesh/face mesh surface count once file post/testmassflow.txt inside_out
#fix massflowout all massflow/mesh/face mesh surface count multiple file post/testmassflow.txt inside_out
# insertion
variable d equal 0.02
create_atoms 1 single -0.55 -0.55 1.1
create_atoms 1 single -0.65 -0.65 1.1
create_atoms 1 single -0.45 -0.45 1.1
create_atoms 1 single -0.45 -0.65 1.1
create_atoms 1 single -0.65 -0.45 1.1
set group all density 3000 diameter $d
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 1000 post/massflow*.vtk id type radius x y z vx vy vz fx fy fz
# total amount of atoms passed through each face
variable totalatomsout1 equal f_massflowout[1][2]
variable totalatomsout2 equal f_massflowout[2][2]
variable totalatomsout3 equal f_massflowout[3][2]
variable totalatomsout4 equal f_massflowout[4][2]
variable totalatomsout5 equal f_massflowout[5][2]
variable totalatomsout6 equal f_massflowout[6][2]
variable totalatomsout7 equal f_massflowout[7][2]
variable totalatomsout8 equal f_massflowout[8][2]
variable totalatomsout9 equal f_massflowout[9][2]
variable totalatomsout10 equal f_massflowout[10][2]
variable totalatomsout11 equal f_massflowout[11][2]
variable totalatomsout12 equal f_massflowout[12][2]
variable totalatomsout13 equal f_massflowout[13][2]
variable totalatomsout14 equal f_massflowout[14][2]
variable totalatomsout15 equal f_massflowout[15][2]
variable totalatomsout16 equal f_massflowout[16][2]
variable totalatomsout17 equal f_massflowout[17][2]
variable totalatomsout18 equal f_massflowout[18][2]
variable totalatomsout19 equal f_massflowout[19][2]
variable totalatomsout20 equal f_massflowout[20][2]
variable totalatomsout21 equal f_massflowout[21][2]
variable totalatomsout22 equal f_massflowout[22][2]
variable totalatomsout23 equal f_massflowout[23][2]
variable totalatomsout24 equal f_massflowout[24][2]
# current amount of atoms passed through each face
variable atomsout1 equal f_massflowout[1][7]
variable atomsout2 equal f_massflowout[2][7]
variable atomsout3 equal f_massflowout[3][7]
variable atomsout4 equal f_massflowout[4][7]
variable atomsout5 equal f_massflowout[5][7]
variable atomsout6 equal f_massflowout[6][7]
variable atomsout7 equal f_massflowout[7][7]
variable atomsout8 equal f_massflowout[8][7]
variable atomsout9 equal f_massflowout[9][7]
variable atomsout10 equal f_massflowout[10][7]
variable atomsout11 equal f_massflowout[11][7]
variable atomsout12 equal f_massflowout[12][7]
variable atomsout13 equal f_massflowout[13][7]
variable atomsout14 equal f_massflowout[14][7]
variable atomsout15 equal f_massflowout[15][7]
variable atomsout16 equal f_massflowout[16][7]
variable atomsout17 equal f_massflowout[17][7]
variable atomsout18 equal f_massflowout[18][7]
variable atomsout19 equal f_massflowout[19][7]
variable atomsout20 equal f_massflowout[20][7]
variable atomsout21 equal f_massflowout[21][7]
variable atomsout22 equal f_massflowout[22][7]
variable atomsout23 equal f_massflowout[23][7]
variable atomsout24 equal f_massflowout[24][7]
variable vxout1 equal f_massflowout[1][8]
variable vxout2 equal f_massflowout[2][8]
variable vxout3 equal f_massflowout[3][8]
variable vxout4 equal f_massflowout[4][8]
variable vxout5 equal f_massflowout[5][8]
variable vxout6 equal f_massflowout[6][8]
variable vxout7 equal f_massflowout[7][8]
variable vxout8 equal f_massflowout[8][8]
variable vxout9 equal f_massflowout[9][8]
variable vxout10 equal f_massflowout[10][8]
variable vxout11 equal f_massflowout[11][8]
variable vxout12 equal f_massflowout[12][8]
variable vxout13 equal f_massflowout[13][8]
variable vxout14 equal f_massflowout[14][8]
variable vxout15 equal f_massflowout[15][8]
variable vxout16 equal f_massflowout[16][8]
variable vxout17 equal f_massflowout[17][8]
variable vxout18 equal f_massflowout[18][8]
variable vxout19 equal f_massflowout[19][8]
variable vxout20 equal f_massflowout[20][8]
variable vxout21 equal f_massflowout[21][8]
variable vxout22 equal f_massflowout[22][8]
variable vxout23 equal f_massflowout[23][8]
variable vxout24 equal f_massflowout[24][8]
variable vyout1 equal f_massflowout[1][9]
variable vyout2 equal f_massflowout[2][9]
variable vyout3 equal f_massflowout[3][9]
variable vyout4 equal f_massflowout[4][9]
variable vyout5 equal f_massflowout[5][9]
variable vyout6 equal f_massflowout[6][9]
variable vyout7 equal f_massflowout[7][9]
variable vyout8 equal f_massflowout[8][9]
variable vyout9 equal f_massflowout[9][9]
variable vyout10 equal f_massflowout[10][9]
variable vyout11 equal f_massflowout[11][9]
variable vyout12 equal f_massflowout[12][9]
variable vyout13 equal f_massflowout[13][9]
variable vyout14 equal f_massflowout[14][9]
variable vyout15 equal f_massflowout[15][9]
variable vyout16 equal f_massflowout[16][9]
variable vyout17 equal f_massflowout[17][9]
variable vyout18 equal f_massflowout[18][9]
variable vyout19 equal f_massflowout[19][9]
variable vyout20 equal f_massflowout[20][9]
variable vyout21 equal f_massflowout[21][9]
variable vyout22 equal f_massflowout[22][9]
variable vyout23 equal f_massflowout[23][9]
variable vyout24 equal f_massflowout[24][9]
variable vzout1 equal f_massflowout[1][10]
variable vzout2 equal f_massflowout[2][10]
variable vzout3 equal f_massflowout[3][10]
variable vzout4 equal f_massflowout[4][10]
variable vzout5 equal f_massflowout[5][10]
variable vzout6 equal f_massflowout[6][10]
variable vzout7 equal f_massflowout[7][10]
variable vzout8 equal f_massflowout[8][10]
variable vzout9 equal f_massflowout[9][10]
variable vzout10 equal f_massflowout[10][10]
variable vzout11 equal f_massflowout[11][10]
variable vzout12 equal f_massflowout[12][10]
variable vzout13 equal f_massflowout[13][10]
variable vzout14 equal f_massflowout[14][10]
variable vzout15 equal f_massflowout[15][10]
variable vzout16 equal f_massflowout[16][10]
variable vzout17 equal f_massflowout[17][10]
variable vzout18 equal f_massflowout[18][10]
variable vzout19 equal f_massflowout[19][10]
variable vzout20 equal f_massflowout[20][10]
variable vzout21 equal f_massflowout[21][10]
variable vzout22 equal f_massflowout[22][10]
variable vzout23 equal f_massflowout[23][10]
variable vzout24 equal f_massflowout[24][10]
variable mstep equal step
fix extratotalatoms all print 1000 "${mstep} ${totalatomsout1} ${totalatomsout2} ${totalatomsout3} ${totalatomsout4} ${totalatomsout5} ${totalatomsout6} ${totalatomsout7} ${totalatomsout8} ${totalatomsout9} ${totalatomsout10} ${totalatomsout11} ${totalatomsout12} ${totalatomsout13} ${totalatomsout14} ${totalatomsout15} ${totalatomsout16} ${totalatomsout17} ${totalatomsout18} ${totalatomsout19} ${totalatomsout20} ${totalatomsout21} ${totalatomsout22} ${totalatomsout23} ${totalatomsout24} " screen no title "step totalatomsout" file post/totalatomsout.txt
fix extraatoms all print 1000 "${mstep} ${atomsout1} ${atomsout2} ${atomsout3} ${atomsout4} ${atomsout5} ${atomsout6} ${atomsout7} ${atomsout8} ${atomsout9} ${atomsout10} ${atomsout11} ${atomsout12} ${atomsout13} ${atomsout14} ${atomsout15} ${atomsout16} ${atomsout17} ${atomsout18} ${atomsout19} ${atomsout20} ${atomsout21} ${atomsout22} ${atomsout23} ${atomsout24} " screen no title "step atomsout" file post/atomsout.txt
fix extravx all print 1000 "${mstep} ${vxout1} ${vxout2} ${vxout3} ${vxout4} ${vxout5} ${vxout6} ${vxout7} ${vxout8} ${vxout9} ${vxout10} ${vxout11} ${vxout12} ${vxout13} ${vxout14} ${vxout15} ${vxout16} ${vxout17} ${vxout18} ${vxout19} ${vxout20} ${vxout21} ${vxout22} ${vxout23} ${vxout24} " screen no title "step vxout" file post/vxout.txt
fix extravy all print 1000 "${mstep} ${vyout1} ${vyout2} ${vyout3} ${vyout4} ${vyout5} ${vyout6} ${vyout7} ${vyout8} ${vyout9} ${vyout10} ${vyout11} ${vyout12} ${vyout13} ${vyout14} ${vyout15} ${vyout16} ${vyout17} ${vyout18} ${vyout19} ${vxout20} ${vyout21} ${vyout22} ${vyout23} ${vyout24} " screen no title "step vyout" file post/vyout.txt
fix extravz all print 1000 "${mstep} ${vzout1} ${vzout2} ${vzout3} ${vzout4} ${vzout5} ${vzout6} ${vzout7} ${vzout8} ${vzout9} ${vzout10} ${vzout11} ${vzout12} ${vzout13} ${vzout14} ${vzout15} ${vzout16} ${vzout17} ${vzout18} ${vzout19} ${vzout20} ${vzout21} ${vzout22} ${vzout23} ${vzout24} " screen no title "step vzout" file post/vzout.txt
run 30000
velocity all zero linear
unfix gravi
fix gravi all gravity 9.81 vector -1.0 -1.0 1.0
run 30000
velocity all zero linear
unfix gravi
fix gravi all gravity 9.81 vector 1.0 1.0 -1.0
run 30000
dump e_data all mesh/vtk 1 post/surface*.vtk face_id
run 1

View File

@ -0,0 +1,380 @@
# measure mass flow through mesh surface (broken down to individual faces)
# additionally measure average temperature and relative layer radii per face
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.2 1.2 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
# material properties required for granular pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.8
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
# pair style
pair_style gran model hertz tangential history
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 1.0 1.0 -1.0
extract_surface meshes/test.vtk file meshes/testsurface.vtk meshes/testextrusion.vtk extrude_length 0.1 min_rad 0.025
fix surface all mesh/surface file meshes/testsurface.vtk type 1 verbose yes cell_data yes
#fix massflowout all massflow/mesh/face mesh surface count once file post/testmassflow.txt inside_out temperature yes chemistry yes
fix massflowout all massflow/mesh/face mesh surface count multiple file post/testmassflow.txt inside_out temperature yes chemistry yes
################################################################################
# heat transfer
fix ftco all property/global thermalConductivity peratomtype 100.
fix ftca all property/global thermalCapacity peratomtype 10.
fix heattransfer all heat/gran initial_temperature 300.
################################################################################
# chemistry
# input fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix partTemp all property/atom partTemp scalar yes no no 1223.15 # Kelvin
fix partNu all property/atom partNu scalar yes no no 0.00004125
fix partRe all property/atom partRe scalar yes no no 2.7500
fix partP all property/atom partP scalar yes no no 101325
# molar fractions
fix X_CO all property/atom X_CO scalar yes no no 0.4
fix X_CO2 all property/atom X_CO2 scalar yes no no 0.0
fix X_N2 all property/atom X_N2 scalar yes no no 0.6
# molecular diffusion coefficient
fix CO_diffCoeff all property/atom CO_diffCoeff scalar yes no no 0.00023
# output fixes for fix chem/shrink/core that would otherwise be provided by fix
# couple/cfd/chemistry
fix Modified_CO all property/atom Modified_CO scalar yes yes no 0.0
fix Modified_CO2 all property/atom Modified_CO2 scalar yes yes no 0.0
fix reactionHeat all property/atom reactionHeat scalar yes no no 0.0
################################################################################
# activate 3-layer unreacted shrinking core model
fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 &
speciesC CO2 molMassC 0.04401 &
scale_reduction_rate 2.0 nevery 25 #screen yes
# chemical properties for unreacted shrinking core model
# w->Fe m->w h->m
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.61 0.31 0.16 0.15
fix tortuosity all property/global tortuosity_all scalar 3
fix pore_diameter all property/global pore_diameter_all scalar 5.5e-7
# define layer properties
fix LayerRelRadii all property/atom relRadii vector yes no no 1.0 0.998 0.995 0.98
# define fix for layer densities
fix LayerDensities all property/global density_all vector 7870. 5740. 5170. 5240.
################################################################################
# insertion
variable d equal 0.02
create_atoms 1 single -0.55 -0.55 1.1
create_atoms 1 single -0.65 -0.65 1.1
create_atoms 1 single -0.45 -0.45 1.1
create_atoms 1 single -0.45 -0.65 1.1
create_atoms 1 single -0.65 -0.45 1.1
set group all density 3000 diameter $d
# apply nve integration
fix integr all nve/sphere
# output settings
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic no extra 0
run 1
dump dmp all custom/vtk 1000 post/massflow*.vtk id type radius x y z vx vy vz fx fy fz
# total amount of atoms passed through each face
variable totalatomsout1 equal f_massflowout[1][2]
variable totalatomsout2 equal f_massflowout[2][2]
variable totalatomsout3 equal f_massflowout[3][2]
variable totalatomsout4 equal f_massflowout[4][2]
variable totalatomsout5 equal f_massflowout[5][2]
variable totalatomsout6 equal f_massflowout[6][2]
variable totalatomsout7 equal f_massflowout[7][2]
variable totalatomsout8 equal f_massflowout[8][2]
variable totalatomsout9 equal f_massflowout[9][2]
variable totalatomsout10 equal f_massflowout[10][2]
variable totalatomsout11 equal f_massflowout[11][2]
variable totalatomsout12 equal f_massflowout[12][2]
variable totalatomsout13 equal f_massflowout[13][2]
variable totalatomsout14 equal f_massflowout[14][2]
variable totalatomsout15 equal f_massflowout[15][2]
variable totalatomsout16 equal f_massflowout[16][2]
variable totalatomsout17 equal f_massflowout[17][2]
variable totalatomsout18 equal f_massflowout[18][2]
variable totalatomsout19 equal f_massflowout[19][2]
variable totalatomsout20 equal f_massflowout[20][2]
variable totalatomsout21 equal f_massflowout[21][2]
variable totalatomsout22 equal f_massflowout[22][2]
variable totalatomsout23 equal f_massflowout[23][2]
variable totalatomsout24 equal f_massflowout[24][2]
# current amount of atoms passed through each face
variable atomsout1 equal f_massflowout[1][7]
variable atomsout2 equal f_massflowout[2][7]
variable atomsout3 equal f_massflowout[3][7]
variable atomsout4 equal f_massflowout[4][7]
variable atomsout5 equal f_massflowout[5][7]
variable atomsout6 equal f_massflowout[6][7]
variable atomsout7 equal f_massflowout[7][7]
variable atomsout8 equal f_massflowout[8][7]
variable atomsout9 equal f_massflowout[9][7]
variable atomsout10 equal f_massflowout[10][7]
variable atomsout11 equal f_massflowout[11][7]
variable atomsout12 equal f_massflowout[12][7]
variable atomsout13 equal f_massflowout[13][7]
variable atomsout14 equal f_massflowout[14][7]
variable atomsout15 equal f_massflowout[15][7]
variable atomsout16 equal f_massflowout[16][7]
variable atomsout17 equal f_massflowout[17][7]
variable atomsout18 equal f_massflowout[18][7]
variable atomsout19 equal f_massflowout[19][7]
variable atomsout20 equal f_massflowout[20][7]
variable atomsout21 equal f_massflowout[21][7]
variable atomsout22 equal f_massflowout[22][7]
variable atomsout23 equal f_massflowout[23][7]
variable atomsout24 equal f_massflowout[24][7]
variable vxout1 equal f_massflowout[1][8]
variable vxout2 equal f_massflowout[2][8]
variable vxout3 equal f_massflowout[3][8]
variable vxout4 equal f_massflowout[4][8]
variable vxout5 equal f_massflowout[5][8]
variable vxout6 equal f_massflowout[6][8]
variable vxout7 equal f_massflowout[7][8]
variable vxout8 equal f_massflowout[8][8]
variable vxout9 equal f_massflowout[9][8]
variable vxout10 equal f_massflowout[10][8]
variable vxout11 equal f_massflowout[11][8]
variable vxout12 equal f_massflowout[12][8]
variable vxout13 equal f_massflowout[13][8]
variable vxout14 equal f_massflowout[14][8]
variable vxout15 equal f_massflowout[15][8]
variable vxout16 equal f_massflowout[16][8]
variable vxout17 equal f_massflowout[17][8]
variable vxout18 equal f_massflowout[18][8]
variable vxout19 equal f_massflowout[19][8]
variable vxout20 equal f_massflowout[20][8]
variable vxout21 equal f_massflowout[21][8]
variable vxout22 equal f_massflowout[22][8]
variable vxout23 equal f_massflowout[23][8]
variable vxout24 equal f_massflowout[24][8]
variable vyout1 equal f_massflowout[1][9]
variable vyout2 equal f_massflowout[2][9]
variable vyout3 equal f_massflowout[3][9]
variable vyout4 equal f_massflowout[4][9]
variable vyout5 equal f_massflowout[5][9]
variable vyout6 equal f_massflowout[6][9]
variable vyout7 equal f_massflowout[7][9]
variable vyout8 equal f_massflowout[8][9]
variable vyout9 equal f_massflowout[9][9]
variable vyout10 equal f_massflowout[10][9]
variable vyout11 equal f_massflowout[11][9]
variable vyout12 equal f_massflowout[12][9]
variable vyout13 equal f_massflowout[13][9]
variable vyout14 equal f_massflowout[14][9]
variable vyout15 equal f_massflowout[15][9]
variable vyout16 equal f_massflowout[16][9]
variable vyout17 equal f_massflowout[17][9]
variable vyout18 equal f_massflowout[18][9]
variable vyout19 equal f_massflowout[19][9]
variable vyout20 equal f_massflowout[20][9]
variable vyout21 equal f_massflowout[21][9]
variable vyout22 equal f_massflowout[22][9]
variable vyout23 equal f_massflowout[23][9]
variable vyout24 equal f_massflowout[24][9]
variable vzout1 equal f_massflowout[1][10]
variable vzout2 equal f_massflowout[2][10]
variable vzout3 equal f_massflowout[3][10]
variable vzout4 equal f_massflowout[4][10]
variable vzout5 equal f_massflowout[5][10]
variable vzout6 equal f_massflowout[6][10]
variable vzout7 equal f_massflowout[7][10]
variable vzout8 equal f_massflowout[8][10]
variable vzout9 equal f_massflowout[9][10]
variable vzout10 equal f_massflowout[10][10]
variable vzout11 equal f_massflowout[11][10]
variable vzout12 equal f_massflowout[12][10]
variable vzout13 equal f_massflowout[13][10]
variable vzout14 equal f_massflowout[14][10]
variable vzout15 equal f_massflowout[15][10]
variable vzout16 equal f_massflowout[16][10]
variable vzout17 equal f_massflowout[17][10]
variable vzout18 equal f_massflowout[18][10]
variable vzout19 equal f_massflowout[19][10]
variable vzout20 equal f_massflowout[20][10]
variable vzout21 equal f_massflowout[21][10]
variable vzout22 equal f_massflowout[22][10]
variable vzout23 equal f_massflowout[23][10]
variable vzout24 equal f_massflowout[24][10]
# current temperature of atoms passed through each face
variable tempout1 equal f_massflowout[1][14]
variable tempout2 equal f_massflowout[2][14]
variable tempout3 equal f_massflowout[3][14]
variable tempout4 equal f_massflowout[4][14]
variable tempout5 equal f_massflowout[5][14]
variable tempout6 equal f_massflowout[6][14]
variable tempout7 equal f_massflowout[7][14]
variable tempout8 equal f_massflowout[8][14]
variable tempout9 equal f_massflowout[9][14]
variable tempout10 equal f_massflowout[10][14]
variable tempout11 equal f_massflowout[11][14]
variable tempout12 equal f_massflowout[12][14]
variable tempout13 equal f_massflowout[13][14]
variable tempout14 equal f_massflowout[14][14]
variable tempout15 equal f_massflowout[15][14]
variable tempout16 equal f_massflowout[16][14]
variable tempout17 equal f_massflowout[17][14]
variable tempout18 equal f_massflowout[18][14]
variable tempout19 equal f_massflowout[19][14]
variable tempout20 equal f_massflowout[20][14]
variable tempout21 equal f_massflowout[21][14]
variable tempout22 equal f_massflowout[22][14]
variable tempout23 equal f_massflowout[23][14]
variable tempout24 equal f_massflowout[24][14]
# current relative layer radii of atoms passed through each face
variable relradw1 equal f_massflowout[1][15]
variable relradw2 equal f_massflowout[2][15]
variable relradw3 equal f_massflowout[3][15]
variable relradw4 equal f_massflowout[4][15]
variable relradw5 equal f_massflowout[5][15]
variable relradw6 equal f_massflowout[6][15]
variable relradw7 equal f_massflowout[7][15]
variable relradw8 equal f_massflowout[8][15]
variable relradw9 equal f_massflowout[9][15]
variable relradw10 equal f_massflowout[10][15]
variable relradw11 equal f_massflowout[11][15]
variable relradw12 equal f_massflowout[12][15]
variable relradw13 equal f_massflowout[13][15]
variable relradw14 equal f_massflowout[14][15]
variable relradw15 equal f_massflowout[15][15]
variable relradw16 equal f_massflowout[16][15]
variable relradw17 equal f_massflowout[17][15]
variable relradw18 equal f_massflowout[18][15]
variable relradw19 equal f_massflowout[19][15]
variable relradw20 equal f_massflowout[20][15]
variable relradw21 equal f_massflowout[21][15]
variable relradw22 equal f_massflowout[22][15]
variable relradw23 equal f_massflowout[23][15]
variable relradw24 equal f_massflowout[24][15]
variable relradm1 equal f_massflowout[1][16]
variable relradm2 equal f_massflowout[2][16]
variable relradm3 equal f_massflowout[3][16]
variable relradm4 equal f_massflowout[4][16]
variable relradm5 equal f_massflowout[5][16]
variable relradm6 equal f_massflowout[6][16]
variable relradm7 equal f_massflowout[7][16]
variable relradm8 equal f_massflowout[8][16]
variable relradm9 equal f_massflowout[9][16]
variable relradm10 equal f_massflowout[10][16]
variable relradm11 equal f_massflowout[11][16]
variable relradm12 equal f_massflowout[12][16]
variable relradm13 equal f_massflowout[13][16]
variable relradm14 equal f_massflowout[14][16]
variable relradm15 equal f_massflowout[15][16]
variable relradm16 equal f_massflowout[16][16]
variable relradm17 equal f_massflowout[17][16]
variable relradm18 equal f_massflowout[18][16]
variable relradm19 equal f_massflowout[19][16]
variable relradm20 equal f_massflowout[20][16]
variable relradm21 equal f_massflowout[21][16]
variable relradm22 equal f_massflowout[22][16]
variable relradm23 equal f_massflowout[23][16]
variable relradm24 equal f_massflowout[24][16]
variable relradh1 equal f_massflowout[1][17]
variable relradh2 equal f_massflowout[2][17]
variable relradh3 equal f_massflowout[3][17]
variable relradh4 equal f_massflowout[4][17]
variable relradh5 equal f_massflowout[5][17]
variable relradh6 equal f_massflowout[6][17]
variable relradh7 equal f_massflowout[7][17]
variable relradh8 equal f_massflowout[8][17]
variable relradh9 equal f_massflowout[9][17]
variable relradh10 equal f_massflowout[10][17]
variable relradh11 equal f_massflowout[11][17]
variable relradh12 equal f_massflowout[12][17]
variable relradh13 equal f_massflowout[13][17]
variable relradh14 equal f_massflowout[14][17]
variable relradh15 equal f_massflowout[15][17]
variable relradh16 equal f_massflowout[16][17]
variable relradh17 equal f_massflowout[17][17]
variable relradh18 equal f_massflowout[18][17]
variable relradh19 equal f_massflowout[19][17]
variable relradh20 equal f_massflowout[20][17]
variable relradh21 equal f_massflowout[21][17]
variable relradh22 equal f_massflowout[22][17]
variable relradh23 equal f_massflowout[23][17]
variable relradh24 equal f_massflowout[24][17]
variable mstep equal step
fix extratotalatoms all print 1000 "${mstep} ${totalatomsout1} ${totalatomsout2} ${totalatomsout3} ${totalatomsout4} ${totalatomsout5} ${totalatomsout6} ${totalatomsout7} ${totalatomsout8} ${totalatomsout9} ${totalatomsout10} ${totalatomsout11} ${totalatomsout12} ${totalatomsout13} ${totalatomsout14} ${totalatomsout15} ${totalatomsout16} ${totalatomsout17} ${totalatomsout18} ${totalatomsout19} ${totalatomsout20} ${totalatomsout21} ${totalatomsout22} ${totalatomsout23} ${totalatomsout24} " screen no title "step totalatomsout" file post/totalatomsout.txt
fix extraatoms all print 1000 "${mstep} ${atomsout1} ${atomsout2} ${atomsout3} ${atomsout4} ${atomsout5} ${atomsout6} ${atomsout7} ${atomsout8} ${atomsout9} ${atomsout10} ${atomsout11} ${atomsout12} ${atomsout13} ${atomsout14} ${atomsout15} ${atomsout16} ${atomsout17} ${atomsout18} ${atomsout19} ${atomsout20} ${atomsout21} ${atomsout22} ${atomsout23} ${atomsout24} " screen no title "step atomsout" file post/atomsout.txt
fix extravx all print 1000 "${mstep} ${vxout1} ${vxout2} ${vxout3} ${vxout4} ${vxout5} ${vxout6} ${vxout7} ${vxout8} ${vxout9} ${vxout10} ${vxout11} ${vxout12} ${vxout13} ${vxout14} ${vxout15} ${vxout16} ${vxout17} ${vxout18} ${vxout19} ${vxout20} ${vxout21} ${vxout22} ${vxout23} ${vxout24} " screen no title "step vxout" file post/vxout.txt
fix extravy all print 1000 "${mstep} ${vyout1} ${vyout2} ${vyout3} ${vyout4} ${vyout5} ${vyout6} ${vyout7} ${vyout8} ${vyout9} ${vyout10} ${vyout11} ${vyout12} ${vyout13} ${vyout14} ${vyout15} ${vyout16} ${vyout17} ${vyout18} ${vyout19} ${vxout20} ${vyout21} ${vyout22} ${vyout23} ${vyout24} " screen no title "step vyout" file post/vyout.txt
fix extravz all print 1000 "${mstep} ${vzout1} ${vzout2} ${vzout3} ${vzout4} ${vzout5} ${vzout6} ${vzout7} ${vzout8} ${vzout9} ${vzout10} ${vzout11} ${vzout12} ${vzout13} ${vzout14} ${vzout15} ${vzout16} ${vzout17} ${vzout18} ${vzout19} ${vzout20} ${vzout21} ${vzout22} ${vzout23} ${vzout24} " screen no title "step vzout" file post/vzout.txt
fix extratemp all print 1000 "${mstep} ${tempout1} ${tempout2} ${tempout3} ${tempout4} ${tempout5} ${tempout6} ${tempout7} ${tempout8} ${tempout9} ${tempout10} ${tempout11} ${tempout12} ${tempout13} ${tempout14} ${tempout15} ${tempout16} ${tempout17} ${tempout18} ${tempout19} ${tempout20} ${tempout21} ${tempout22} ${tempout23} ${tempout24} " screen no title "step temp" file post/tempout.txt
fix extrarelradw all print 1000 "${mstep} ${relradw1} ${relradw2} ${relradw3} ${relradw4} ${relradw5} ${relradw6} ${relradw7} ${relradw8} ${relradw9} ${relradw10} ${relradw11} ${relradw12} ${relradw13} ${relradw14} ${relradw15} ${relradw16} ${relradw17} ${relradw18} ${relradw19} ${relradw20} ${relradw21} ${relradw22} ${relradw23} ${relradw24} " screen no title "step relradw" file post/relradw.txt
fix extrarelradm all print 1000 "${mstep} ${relradm1} ${relradm2} ${relradm3} ${relradm4} ${relradm5} ${relradm6} ${relradm7} ${relradm8} ${relradm9} ${relradm10} ${relradm11} ${relradm12} ${relradm13} ${relradm14} ${relradm15} ${relradm16} ${relradm17} ${relradm18} ${relradm19} ${relradm20} ${relradm21} ${relradm22} ${relradm23} ${relradm24} " screen no title "step relradm" file post/relradm.txt
fix extrarelradh all print 1000 "${mstep} ${relradh1} ${relradh2} ${relradh3} ${relradh4} ${relradh5} ${relradh6} ${relradh7} ${relradh8} ${relradh9} ${relradh10} ${relradh11} ${relradh12} ${relradh13} ${relradh14} ${relradh15} ${relradh16} ${relradh17} ${relradh18} ${relradh19} ${relradh20} ${relradh21} ${relradh22} ${relradh23} ${relradh24} " screen no title "step relradh" file post/relradh.txt
run 30000
set group all property/atom Temp 350
velocity all zero linear
unfix gravi
fix gravi all gravity 9.81 vector -1.0 -1.0 1.0
run 30000
set atom 1 property/atom Temp 320
set atom 2 property/atom Temp 330
velocity all zero linear
unfix gravi
fix gravi all gravity 9.81 vector 1.0 1.0 -1.0
run 30000
dump e_data all mesh/vtk 1 post/surface*.vtk face_id
run 1

View File

@ -0,0 +1,50 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 27 float
-0.5 -0.5 0.0
0.0 -0.5 0.0
0.5 -0.5 0.0
-0.5 0.0 0.0
0.0 0.0 0.0
0.5 0.0 0.0
-0.5 0.5 0.0
0.0 0.5 0.0
0.5 0.5 0.0
-0.5 -0.5 0.5
0.0 -0.5 0.5
0.5 -0.5 0.5
-0.5 0.0 0.5
0.0 0.0 0.5
0.5 0.0 0.5
-0.5 0.5 0.5
0.0 0.5 0.5
0.5 0.5 0.5
-0.5 -0.5 1.0
0.0 -0.5 1.0
0.5 -0.5 1.0
-0.5 0.0 1.0
0.0 0.0 1.0
0.5 0.0 1.0
-0.5 0.5 1.0
0.0 0.5 1.0
0.5 0.5 1.0
CELLS 8 72
8 0 1 4 3 9 10 13 12
8 1 2 5 4 10 11 14 13
8 3 4 7 6 12 13 16 15
8 4 5 8 7 13 14 17 16
8 9 10 13 12 18 19 22 21
8 10 11 14 13 19 20 23 22
8 12 13 16 15 21 22 25 24
8 13 14 17 16 22 23 26 25
CELL_TYPES 8
12
12
12
12
12
12
12
12

View File

@ -0,0 +1,14 @@
{
"runs": [
{
"name" : "massflow_mesh_face",
"input_script" : "in.massflow_mesh_face",
"type" : "serial"
},
{
"name" : "massflow_mesh_face_ex",
"input_script": "in.massflow_mesh_face_ex",
"type": "serial"
}
]
}

View File

@ -44,14 +44,14 @@
"type" : "serial",
"data" : {
"series" : [
{"name" : "rhoeff1", "file" : "rho_1.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff1", "file" : "rho_1.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm1", "file" : "Aterm_1.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm1", "file" : "Bterm_1.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm1", "file" : "MassTerm_1.dat", "columns" : ["time", "mt"]},
{"name" : "dmA1", "file" : "dmdot_1.dat", "columns" : ["time", "dmA1", "dmA2", "dmA3"]},
{"name" : "fracRed1", "file" : "fr_d1_1.dat", "columns" : ["time", "fr1", "fr2", "fr3", "frtot"]},
{"name" : "relrad1", "file" : "relRadii_1.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]},
{"name" : "rhoeff2", "file" : "rho_2.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff2", "file" : "rho_2.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm2", "file" : "Aterm_2.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm2", "file" : "Bterm_2.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm2", "file" : "MassTerm_2.dat", "columns" : ["time", "mt"]},
@ -60,8 +60,8 @@
{"name" : "relrad2", "file" : "relRadii_2.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]}
],
"plots" : [
{"name" : "plot_rhoeff_atom1", "title" : "Effective densities (atom 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff","rhoeff1.rhoeffl1","rhoeff1.rhoeffl2","rhoeff1.rhoeffl3","rhoeff1.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_atom2", "title" : "Effective densities (atom 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff","rhoeff2.rhoeffl1","rhoeff2.rhoeffl2","rhoeff2.rhoeffl3","rhoeff2.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_atom1", "title" : "Effective density (atom 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_rhoeff_atom2", "title" : "Effective density (atom 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_aterm_atom1", "title" : "Reaction resistance terms (atom 1)", "xdata" : "aterm1.time", "ydata" : ["aterm1.a1","aterm1.a2","aterm1.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_aterm_atom2", "title" : "Reaction resistance terms (atom 2)", "xdata" : "aterm2.time", "ydata" : ["aterm2.a1","aterm2.a2","aterm2.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_bterm_atom1", "title" : "Diffusion resistance terms (atom 1)", "xdata" : "bterm1.time", "ydata" : ["bterm1.b1","bterm1.b2","bterm1.b3"], "xlabel" : "Time [s]", "ylabel" : "Diffusion resistance", "legend" : ["B_Fe","B_w","B_m"]},
@ -84,14 +84,14 @@
"type" : "serial",
"data" : {
"series" : [
{"name" : "rhoeff1", "file" : "rho_g1.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff1", "file" : "rho_g1.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm1", "file" : "Aterm_g1.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm1", "file" : "Bterm_g1.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm1", "file" : "MassTerm_g1.dat", "columns" : ["time", "mt"]},
{"name" : "dmA1", "file" : "dmdot_g1.dat", "columns" : ["time", "dmA1", "dmA2", "dmA3"]},
{"name" : "fracRed1", "file" : "fr_d1_g1.dat", "columns" : ["time", "fr1", "fr2", "fr3", "frtot"]},
{"name" : "relrad1", "file" : "relRadii_g1.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]},
{"name" : "rhoeff2", "file" : "rho_g2.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff2", "file" : "rho_g2.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm2", "file" : "Aterm_g2.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm2", "file" : "Bterm_g2.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm2", "file" : "MassTerm_g2.dat", "columns" : ["time", "mt"]},
@ -100,8 +100,8 @@
{"name" : "relrad2", "file" : "relRadii_g2.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]}
],
"plots" : [
{"name" : "plot_rhoeff_group1", "title" : "Effective densities (group 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff","rhoeff1.rhoeffl1","rhoeff1.rhoeffl2","rhoeff1.rhoeffl3","rhoeff1.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_group2", "title" : "Effective densities (group 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff","rhoeff2.rhoeffl1","rhoeff2.rhoeffl2","rhoeff2.rhoeffl3","rhoeff2.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_group1", "title" : "Effective density (group 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_rhoeff_group2", "title" : "Effective density (group 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_aterm_group1", "title" : "Reaction resistance terms (group 1)", "xdata" : "aterm1.time", "ydata" : ["aterm1.a1","aterm1.a2","aterm1.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_aterm_group2", "title" : "Reaction resistance terms (group 2)", "xdata" : "aterm2.time", "ydata" : ["aterm2.a1","aterm2.a2","aterm2.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_bterm_group1", "title" : "Diffusion resistance terms (group 1)", "xdata" : "bterm1.time", "ydata" : ["bterm1.b1","bterm1.b2","bterm1.b3"], "xlabel" : "Time [s]", "ylabel" : "Diffusion resistance", "legend" : ["B_Fe","B_w","B_m"]},
@ -124,14 +124,14 @@
"type" : "serial",
"data" : {
"series" : [
{"name" : "rhoeff1", "file" : "rho_g1.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff1", "file" : "rho_g1.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm1", "file" : "Aterm_g1.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm1", "file" : "Bterm_g1.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm1", "file" : "MassTerm_g1.dat", "columns" : ["time", "mt"]},
{"name" : "dmA1", "file" : "dmdot_g1.dat", "columns" : ["time", "dmA1", "dmA2", "dmA3"]},
{"name" : "fracRed1", "file" : "fr_d1_g1.dat", "columns" : ["time", "fr1", "fr2", "fr3", "frtot"]},
{"name" : "relrad1", "file" : "relRadii_g1.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]},
{"name" : "rhoeff2", "file" : "rho_g2.dat", "columns" : ["time", "rhoeff", "rhoeffl1", "rhoeffl2", "rhoeffl3", "rhoeffl4"]},
{"name" : "rhoeff2", "file" : "rho_g2.dat", "columns" : ["time", "rhoeff"]},
{"name" : "aterm2", "file" : "Aterm_g2.dat", "columns" : ["time", "a1", "a2", "a3"]},
{"name" : "bterm2", "file" : "Bterm_g2.dat", "columns" : ["time", "b1", "b2", "b3"]},
{"name" : "mterm2", "file" : "MassTerm_g2.dat", "columns" : ["time", "mt"]},
@ -140,8 +140,8 @@
{"name" : "relrad2", "file" : "relRadii_g2.dat", "columns" : ["time", "rr1", "rr2", "rr3", "rr4"]}
],
"plots" : [
{"name" : "plot_rhoeff_group1", "title" : "Effective densities (group 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff","rhoeff1.rhoeffl1","rhoeff1.rhoeffl2","rhoeff1.rhoeffl3","rhoeff1.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_group2", "title" : "Effective densities (group 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff","rhoeff2.rhoeffl1","rhoeff2.rhoeffl2","rhoeff2.rhoeffl3","rhoeff2.rhoeffl4"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff","rhoeffFe","rhoeffW","rhoeffM","rhoeffH"]},
{"name" : "plot_rhoeff_group1", "title" : "Effective density (group 1)", "xdata" : "rhoeff1.time", "ydata" : ["rhoeff1.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_rhoeff_group2", "title" : "Effective density (group 2)", "xdata" : "rhoeff2.time", "ydata" : ["rhoeff2.rhoeff"], "xlabel" : "Time [s]", "ylabel" : "Effective density [kg/m3]", "legend" : ["rhoeff"]},
{"name" : "plot_aterm_group1", "title" : "Reaction resistance terms (group 1)", "xdata" : "aterm1.time", "ydata" : ["aterm1.a1","aterm1.a2","aterm1.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_aterm_group2", "title" : "Reaction resistance terms (group 2)", "xdata" : "aterm2.time", "ydata" : ["aterm2.a1","aterm2.a2","aterm2.a3"], "xlabel" : "Time [s]", "ylabel" : "Reaction resistance", "legend" : ["A_w->Fe","A_m->w","A_h->m"]},
{"name" : "plot_bterm_group1", "title" : "Diffusion resistance terms (group 1)", "xdata" : "bterm1.time", "ydata" : ["bterm1.b1","bterm1.b2","bterm1.b3"], "xlabel" : "Time [s]", "ylabel" : "Diffusion resistance", "legend" : ["B_Fe","B_w","B_m"]},

View File

@ -0,0 +1,67 @@
atom_style granular
atom_modify map array
boundary m m m
newton off
communicate single vel yes
#processors 2 1 1
units si
region reg block -1 1 -1 1 -0.05 2.85 units box
create_box 1 reg
neighbor 0.01 bin
neigh_modify delay 0
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
fix m5 all property/global characteristicVelocity scalar 2.
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 500
fix customProp all property/atom customProperty scalar yes no no 0.
#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
#particle distributions
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2500 radius constant 0.015
fix pts2 all particletemplate/sphere 1 atom_type 1 density constant 2500 radius constant 0.025
fix pdd1 all particledistribution/discrete 1. 2 pts1 0.3 pts2 0.7
#region and insertion
fix ins_mesh all mesh/surface file meshes/face.stl type 1 scale 0.005
group nve_group region reg
fix ins1 nve_group insert/stream seed 100001 distributiontemplate pdd1 nparticles 2500 &
vel constant 0. -0.5 -2. particlerate 500 set_property customProperty 1. &
overlapcheck yes insertion_face ins_mesh extrude_length 0.6
fix ins2 nve_group insert/stream seed 100001 distributiontemplate pdd1 nparticles 2500 &
vel constant 0. -0.5 -2. particlerate 500 set_property customProperty 2. &
overlapcheck yes insertion_face ins_mesh extrude_length 0.6
#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere
#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
dump dmp all custom 800 post/dump*.stream_property id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius f_customProp
#insert particles
run 100000

View File

@ -0,0 +1,68 @@
atom_style granular
atom_modify map array
boundary m m m
newton off
communicate single vel yes
#processors 2 1 1
units si
region reg block -1 1 -1 1 -0.05 2.85 units box
create_box 1 reg
neighbor 0.01 bin
neigh_modify delay 0
#Material properties required for new pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
fix m5 all property/global characteristicVelocity scalar 2.
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 500
fix customProp all property/atom/lammps i_customProperty
#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
#particle distributions
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2500 radius constant 0.015
fix pts2 all particletemplate/sphere 1 atom_type 1 density constant 2500 radius constant 0.025
fix pdd1 all particledistribution/discrete 1. 2 pts1 0.3 pts2 0.7
#region and insertion
fix ins_mesh all mesh/surface file meshes/face.stl type 1 scale 0.005
group nve_group region reg
fix ins1 nve_group insert/stream seed 100001 distributiontemplate pdd1 nparticles 2500 &
vel constant 0. -0.5 -2. particlerate 500 set_property i_customProperty 1 &
overlapcheck yes insertion_face ins_mesh extrude_length 0.6
fix ins2 nve_group insert/stream seed 100001 distributiontemplate pdd1 nparticles 2500 &
vel constant 0. -0.5 -2. particlerate 500 set_property i_customProperty 2 &
overlapcheck yes insertion_face ins_mesh extrude_length 0.6
#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere
#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
compute custom all property/atom i_customProperty
dump dmp all custom 800 post/dump*.stream_property_lammps id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius c_custom
#insert particles
run 100000

View File

@ -4,6 +4,16 @@
"name" : "serial",
"input_script" : "in.insert_stream",
"type" : "serial"
},
{
"name" : "stream_property",
"input_script" : "in.insert_stream_property",
"type" : "serial"
},
{
"name" : "stream_property_lammps",
"input_script" : "in.insert_stream_property_lammps",
"type" : "serial"
}
]
}

View File

@ -0,0 +1,61 @@
# vtk DataFile Version 3.0
LIGGGHTS grid/mesh
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 27 float
-0.5 -0.5 0.0
0.0 -0.5 0.0
0.5 -0.5 0.0
-0.5 0.0 0.0
0.0 0.0 0.0
0.5 0.0 0.0
-0.5 0.5 0.0
0.0 0.5 0.0
0.5 0.5 0.0
-0.5 -0.5 0.5
0.0 -0.5 0.5
0.5 -0.5 0.5
-0.5 0.0 0.5
0.0 0.0 0.5
0.5 0.0 0.5
-0.5 0.5 0.5
0.0 0.5 0.5
0.5 0.5 0.5
-0.5 -0.5 1.0
0.0 -0.5 1.0
0.5 -0.5 1.0
-0.5 0.0 1.0
0.0 0.0 1.0
0.5 0.0 1.0
-0.5 0.5 1.0
0.0 0.5 1.0
0.5 0.5 1.0
CELLS 8 72
8 0 1 4 3 9 10 13 12
8 1 2 5 4 10 11 14 13
8 3 4 7 6 12 13 16 15
8 4 5 8 7 13 14 17 16
8 9 10 13 12 18 19 22 21
8 10 11 14 13 19 20 23 22
8 12 13 16 15 21 22 25 24
8 13 14 17 16 22 23 26 25
CELL_TYPES 8
12
12
12
12
12
12
12
12
CELL_DATA 8
SCALARS cell_id int 1
LOOKUP_TABLE default
0
1
2
3
4
5
6
7

View File

@ -0,0 +1,59 @@
# hex mesh region example
atom_style granular
atom_modify map array
boundary m m m
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.05 1.05 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
#Material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
#pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
#region and insertion
variable d equal 0.05
region mesh mesh/hex file data/test.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. units box
group nve_group region reg
lattice sc $d
create_atoms 1 region mesh
set group all density 2000 diameter $d
#apply nve integration
fix integr nve_group nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 800 post/hexmesh*.vtk id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
run 1000

View File

@ -0,0 +1,63 @@
# hex mesh region example
atom_style granular
atom_modify map array
boundary m m m
newton off
communicate single vel yes
units si
region reg block -1 1 -1 1 -0.05 1.05 units box
create_box 1 reg
neighbor 0.002 bin
neigh_modify delay 0
#material properties required for pair styles
fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
fix m5 all property/global characteristicVelocity scalar 2.
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 500
#pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *
timestep 0.00001
fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
#region and insertion
variable r equal 0.025
region mesh mesh/hex file data/test.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. units box
group nve_group region reg
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2000 radius constant ${r}
fix pdd1 all particledistribution/discrete 5331 1 pts1 1.0
fix ins all insert/pack seed 100001 distributiontemplate pdd1 vel constant 0. 0. 0.0 &
insert_every once overlapcheck yes all_in no volumefraction_region 0.3 region mesh
#apply nve integration
fix integr nve_group nve/sphere
#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
#insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 800 post/hexmeshrand*.vtk id type x y z vx vy vz fx fy fz omegax omegay omegaz radius
run 1000

View File

@ -0,0 +1,9 @@
{
"runs": [
{
"name" : "serial",
"input_script" : "in.mesh_hex",
"type" : "serial"
}
]
}

View File

@ -0,0 +1 @@
liggghts < in.mesh_hex

View File

@ -5,15 +5,15 @@ OPTION(USE_SUPERQUADRIC "Superquadric particles" OFF)
OPTION(USE_OPENMP "OpenMP parallelization" OFF)
OPTION(TESTING "TESTING" OFF)
SET(LIGGGHTS_MAJOR_VERSION 20)
SET(LIGGGHTS_MINOR_VERSION 09)
SET(LIGGGHTS_MAJOR_VERSION 21)
SET(LIGGGHTS_MINOR_VERSION 03)
SET(LIGGGHTS_PATCH_VERSION 0)
SET(LIGGGHTS_VERSION ${LIGGGHTS_MAJOR_VERSION}.${LIGGGHTS_MINOR_VERSION}.${LIGGGHTS_PATCH_VERSION})
MESSAGE(STATUS "${LIGGGHTS_VERSION}")
IF (NOT CMAKE_CXX_FLAGS)
IF(CMAKE_COMPILER_IS_GNUCXX)
SET(CMAKE_CXX_FLAGS "-O2 -funroll-loops -fstrict-aliasing -Wall -Wextra -Wno-unused-result -Wno-unused-parameter -Wno-literal-suffix -std=c++11")
SET(CMAKE_CXX_FLAGS "-O2 -funroll-loops -fstrict-aliasing -Wall -Wextra -Wno-unused-result -Wno-unused-parameter -Wno-literal-suffix -Wno-cast-function-type -std=c++11")
ENDIF()
ENDIF()

View File

@ -8,7 +8,7 @@ SHELL = /bin/sh
CC = mpic++
CCFLAGS = -O2 \
-funroll-loops -fstrict-aliasing -Wall -Wno-unused-result -Wno-literal-suffix
-funroll-loops -fstrict-aliasing -Wall -Wno-unused-result -Wno-literal-suffix -Wno-cast-function-type
SHFLAGS =
DEPFLAGS = -M

View File

@ -67,7 +67,7 @@ JPG_LIB =
VTK_INC = -I/usr/include/vtk-5.8
VTK_PATH =
VTK_LIB = -lvtkCommon -lvtkFiltering -lvtkIO
VTK_LIB = -lvtkCommon -lvtkFiltering -lvtkGraphics -lvtkIO
# ---------------------------------------------------------------------
# build rules and dependencies

View File

@ -67,7 +67,7 @@ JPG_LIB =
VTK_INC = -I/usr/include/vtk-6.0
VTK_PATH =
VTK_LIB = -lvtkCommonCore-6.0 -lvtkIOCore-6.0 -lvtkIOXML-6.0 -lvtkIOLegacy-6.0 -lvtkCommonDataModel-6.0
VTK_LIB = -lvtkCommonCore-6.0 -lvtkIOCore-6.0 -lvtkIOXML-6.0 -lvtkIOLegacy-6.0 -lvtkCommonDataModel-6.0 -lvtkFiltersCore-6.0 -lvtkFiltersVerdict-6.0 -lvtkFiltersGeometry-6.0 -lvtkCommonExecutionModel-6.0
# ---------------------------------------------------------------------
# build rules and dependencies

View File

@ -0,0 +1,112 @@
# openmpi = Ubuntu 16.04, mpic++, OpenMPI-1.10
SHELL = /bin/sh
# ---------------------------------------------------------------------
# compiler/linker settings
# specify flags and libraries needed for your compiler
CC = mpic++
CCFLAGS = -O2 \
-std=c++11 -funroll-loops -fstrict-aliasing -Wall -Wno-unused-result -Wno-literal-suffix
DEPFLAGS = -M
LINK = mpic++
LINKFLAGS = -O2
LIB = -lstdc++
ARCHIVE = ar
ARFLAGS = -rcsv
SIZE = size
# ---------------------------------------------------------------------
# LAMMPS-specific settings
# specify settings for LAMMPS features you will use
# if you change any -D setting, do full re-compile after "make clean"
# LAMMPS ifdef settings, OPTIONAL
# see possible settings in doc/Section_start.html#2_2 (step 4)
LMP_INC = -DLAMMPS_GZIP -DLAMMPS_VTK
# MPI library, REQUIRED
# see discussion in doc/Section_start.html#2_2 (step 5)
# can point to dummy MPI library in src/STUBS as in Makefile.serial
# INC = path for mpi.h, MPI compiler settings
# PATH = path for MPI library
# LIB = name of MPI library
MPI_INC =
MPI_PATH =
MPI_LIB =
# FFT library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 6)
# can be left blank to use provided KISS FFT library
# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings
# PATH = path for FFT library
# LIB = name of FFT library
FFT_INC =
FFT_PATH =
FFT_LIB =
# JPEG and/or PNG library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 7)
# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC
# INC = path(s) for jpeglib.h and/or png.h
# PATH = path(s) for JPEG library and/or PNG library
# LIB = name(s) of JPEG library and/or PNG library
JPG_INC =
JPG_PATH =
JPG_LIB =
# VTK library, OPTIONAL
# INC = path for VTK header files
# PATH = path for VTK library
# LIB = name of VTK library
VTK_INC = -I/usr/include/vtk-6.2
VTK_PATH =
VTK_LIB = -lvtkCommonCore-6.2 -lvtkCommonDataModel-6.2 -lvtkCommonExecutionModel-6.2 \
-lvtkIOCore-6.2 -lvtkIOXML-6.2 -lvtkIOParallelXML-6.2 -lvtkIOLegacy-6.2 \
-lvtkFiltersCore-6.2 -lvtkFiltersVerdict-6.2 -lvtkFiltersGeometry-6.2
# ---------------------------------------------------------------------
# build rules and dependencies
# no need to edit this section
include Makefile.package.settings
include Makefile.package
EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(VTK_INC) $(PKG_SYSINC)
EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(VTK_PATH) $(PKG_SYSPATH)
EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(VTK_LIB) $(PKG_SYSLIB)
# Path to src files
vpath %.cpp ..
vpath %.h ..
# Link target
$(EXE): $(OBJ)
$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
$(SIZE) $(EXE)
# Library target
lib: $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
# Compilation rules
%.o:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) -c $<
%.d:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
# Individual dependencies
DEPENDS = $(OBJ:.o=.d)
sinclude $(DEPENDS)

View File

@ -0,0 +1,113 @@
# openmpi = Ubuntu 18.04, mpic++, OpenMPI-1.10
SHELL = /bin/sh
# ---------------------------------------------------------------------
# compiler/linker settings
# specify flags and libraries needed for your compiler
CC = mpic++
CCFLAGS = -O2 \
-std=c++11 -funroll-loops -fstrict-aliasing -Wall -Wno-unused-result -Wno-literal-suffix -Wno-cast-function-type
DEPFLAGS = -M
LINK = mpic++
LINKFLAGS = -O2
LIB = -lstdc++
ARCHIVE = ar
ARFLAGS = -rcsv
SIZE = size
# ---------------------------------------------------------------------
# LAMMPS-specific settings
# specify settings for LAMMPS features you will use
# if you change any -D setting, do full re-compile after "make clean"
# LAMMPS ifdef settings, OPTIONAL
# see possible settings in doc/Section_start.html#2_2 (step 4)
LMP_INC = -DLAMMPS_GZIP -DLAMMPS_VTK
# MPI library, REQUIRED
# see discussion in doc/Section_start.html#2_2 (step 5)
# can point to dummy MPI library in src/STUBS as in Makefile.serial
# INC = path for mpi.h, MPI compiler settings
# PATH = path for MPI library
# LIB = name of MPI library
MPI_INC =
MPI_PATH =
MPI_LIB =
# FFT library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 6)
# can be left blank to use provided KISS FFT library
# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings
# PATH = path for FFT library
# LIB = name of FFT library
FFT_INC =
FFT_PATH =
FFT_LIB =
# JPEG and/or PNG library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 7)
# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC
# INC = path(s) for jpeglib.h and/or png.h
# PATH = path(s) for JPEG library and/or PNG library
# LIB = name(s) of JPEG library and/or PNG library
JPG_INC =
JPG_PATH =
JPG_LIB =
# VTK library, OPTIONAL
# INC = path for VTK header files
# PATH = path for VTK library
# LIB = name of VTK library
VTK_INC = -I/usr/include/vtk-6.3
VTK_PATH =
VTK_LIB = -lvtkCommonCore-6.3 -lvtkCommonDataModel-6.3 -lvtkCommonExecutionModel-6.3 \
-lvtkIOCore-6.3 -lvtkIOXML-6.3 -lvtkIOParallelXML-6.3 -lvtkIOLegacy-6.3 \
-lvtkIOPLY-6.3 -lvtkIOGeometry-6.3 -lvtkFiltersCore-6.3 -lvtkFiltersVerdict-6.3 \
-lvtkFiltersGeometry-6.3 -lvtkFiltersModeling-6.3
# ---------------------------------------------------------------------
# build rules and dependencies
# no need to edit this section
include Makefile.package.settings
include Makefile.package
EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(VTK_INC) $(PKG_SYSINC)
EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(VTK_PATH) $(PKG_SYSPATH)
EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(VTK_LIB) $(PKG_SYSLIB)
# Path to src files
vpath %.cpp ..
vpath %.h ..
# Link target
$(EXE): $(OBJ)
$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
$(SIZE) $(EXE)
# Library target
lib: $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
# Compilation rules
%.o:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) -c $<
%.d:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
# Individual dependencies
DEPENDS = $(OBJ:.o=.d)
sinclude $(DEPENDS)

View File

@ -55,21 +55,21 @@ namespace LAMMPS_NS
virtual void scale(double factor) = 0;
// linear move w/ total and incremental displacement
virtual void move(double *vecTotal, double *vecIncremental) = 0;
virtual void move(const double *vecTotal, const double *vecIncremental) = 0;
// linear move w/ incremental displacement
virtual void move(double *vecIncremental) = 0;
virtual void move(const double *vecIncremental) = 0;
// rotation w/ total and incremental displacement
// calls rotate(double *totalQuat,double *dQuat,double *displacement)
virtual void rotate(double totalAngle, double dAngle, double *axis, double *p) = 0;
virtual void rotate(double totalAngle, double dAngle, const double *axis, const double *p) = 0;
// rotation w/ incremental displacement
// calls rotate(double *dQuat,double *displacement)
virtual void rotate(double dAngle, double *axis, double *p) = 0;
virtual void rotate(double dAngle, const double *axis, const double *p) = 0;
// rotation using quaternions
virtual void rotate(double *totalQ, double *dQ,double *origin) = 0;
virtual void rotate(const double *totalQ, const double *dQ, const double *origin) = 0;
// initialize movement
virtual bool registerMove(bool _scale, bool _translate, bool _rotate) = 0;

View File

@ -7,7 +7,8 @@
Christoph Kloss, christoph.kloss@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright 2012-2015 DCS Computing GmbH, Linz
Copyright 2015- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
@ -23,12 +24,14 @@
Contributing authors:
Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
Philippe Seil (JKU Linz)
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_ASSOCIATIVE_POINTER_ARRAY_H
#define LMP_ASSOCIATIVE_POINTER_ARRAY_H
#include <string.h>
#include <string>
#include <map>
#include "memory.h"
namespace LAMMPS_NS
@ -77,10 +80,11 @@ class AssociativePointerArray
inline void storeOrig(const char *_id,class AssociativePointerArray &orig);
inline bool reset(class AssociativePointerArray &orig);
inline bool reset(const char *_id,class AssociativePointerArray &orig);
inline void setToDefault(int n);
void rotate(double *dQ);
void move(double *delta);
void moveElement(int i,double *delta);
void rotate(const double *dQ);
void move(const double *delta);
void moveElement(int i,const double *delta);
void scale(double factor);
inline int bufSize(int operation,bool scale,bool translate,bool rotate) const;
@ -97,15 +101,13 @@ class AssociativePointerArray
inline int pushElemToBuffer(int n, double *buf, int operation,bool scale,bool translate, bool rotate);
inline int popElemFromBuffer(double *buf, int operation,bool scale,bool translate, bool rotate);
int idToIndex(const char *_id);
void indexToId(int index, char *_id);
bool hasId(const char *_id);
private:
T **content_;
int numElem_, maxElem_;
void growArrays();
std::map<std::string, T*> content_;
typedef typename std::map<std::string, T*>:: iterator content_iterator;
typedef typename std::map<std::string, T*>:: const_iterator content_const_iterator;
};
// *************************************

View File

@ -7,7 +7,8 @@
Christoph Kloss, christoph.kloss@cfdem.com
Copyright 2009-2012 JKU Linz
Copyright 2012- DCS Computing GmbH, Linz
Copyright 2012-2015 DCS Computing GmbH, Linz
Copyright 2015- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
@ -24,6 +25,7 @@
Christoph Kloss (JKU Linz, DCS Computing GmbH, Linz)
Philippe Seil (JKU Linz)
Richard Berger (JKU Linz)
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_ASSOCIATIVE_POINTER_ARRAY_I_H
@ -35,19 +37,14 @@
template<typename T>
AssociativePointerArray<T>::AssociativePointerArray()
: content_(0), numElem_(0), maxElem_(1)
{
content_ = new T*[1];
content_[0] = 0;
}
template<typename T>
AssociativePointerArray<T>::~AssociativePointerArray()
{
for(int i = 0; i < numElem_; i++)
delete content_[i];
delete[] content_;
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
delete it->second;
}
/* ----------------------------------------------------------------------
@ -57,17 +54,9 @@
template<typename T> template<typename U>
U* AssociativePointerArray<T>::add(const char *_id, const char* _comm, const char* _ref, const char *_restart, int _scalePower)
{
if(numElem_ == maxElem_)
growArrays();
content_[_id] = static_cast<T*>(new U(_id,_comm,_ref,_restart,_scalePower));
content_[numElem_] = static_cast<T*>(new U(_id,_comm,_ref,_restart,_scalePower));
numElem_++;
/*NP
printf("numElem_ %d\n",numElem_);
for(int i=0;i<numElem_;i++)
printf(" %d %s\n",i,content_[i]->id_);
*/
return static_cast<U*>(content_[numElem_-1]);
return static_cast<U*>(content_[_id]);
}
/* ----------------------------------------------------------------------
@ -77,15 +66,12 @@
template<typename T>
void AssociativePointerArray<T>::remove(const char *_id)
{
int index = idToIndex(_id);
if(index == -1) return;
numElem_--;
delete content_[index];
if(numElem_ > 0)
content_[index] = content_[numElem_];
content_iterator it = content_.find(_id);
if(it != content_.end())
{
delete it->second;
content_.erase(it);
}
}
/* ----------------------------------------------------------------------
@ -95,9 +81,9 @@
template<typename T>
bool AssociativePointerArray<T>::sameLength(int _len)
{
for(int i = 0; i < numElem_; i++)
if(content_[i]->size() != _len)
return false;
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
if(it->second->size() != _len)
return false;
return true;
}
@ -108,9 +94,9 @@
template<typename T> template<typename U>
U* AssociativePointerArray<T>::getPointerById(const char *_id)
{
int ind = idToIndex(_id);
if(ind != -1)
return getPointerByIndex<U>(ind);
content_iterator it = content_.find(_id);
if(it != content_.end())
return dynamic_cast<U*>(it->second);
else
return 0;
}
@ -118,9 +104,9 @@
template<typename T>
T* AssociativePointerArray<T>::getBasePointerById(const char *_id)
{
int ind = idToIndex(_id);
if(ind != -1)
return getBasePointerByIndex(ind);
content_iterator it = content_.find(_id);
if(it != content_.end())
return it->second;
else
return 0;
}
@ -143,41 +129,22 @@
memory management
------------------------------------------------------------------------- */
template<typename T>
void AssociativePointerArray<T>::growArrays()
{
// for(int i=0;i<numElem_+1;i++)
// printf("%d %s %d\n",i,id_[i], strcmp(id_[i],"v"));
T ** tmp = new T*[++maxElem_];
for(int i = 0; i < numElem_; i++)
tmp[i] = content_[i];
delete[] content_;
content_ = tmp;
//for(int i=0;i<numElem_+1;i++)
// printf("%d %s %d\n",i,id_[i], strcmp(id_[i],"v"));
}
template<typename T>
void AssociativePointerArray<T>::grow(int to)
{
int by;
for(int i = 0; i < maxElem_; i++)
{
by = to - getBasePointerByIndex(i)->size();
if(by > 0)
getBasePointerByIndex(i)->addUninitialized(by);
}
{
int by;
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
{
by = to - it->second->size();
if(by > 0)
it->second->addUninitialized(by);
}
}
template<typename T>
int AssociativePointerArray<T>::size() const
{
return numElem_;
return content_.size();
}
/* ----------------------------------------------------------------------
@ -187,8 +154,8 @@
template<typename T>
void AssociativePointerArray<T>::copyElement(int from, int to)
{
for(int i=0;i<numElem_;i++)
content_[i]->copy(from,to);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->copy(from,to);
}
/* ----------------------------------------------------------------------
@ -198,8 +165,8 @@
template<typename T>
void AssociativePointerArray<T>::addUninitializedElement()
{
for(int i=0;i<numElem_;i++)
content_[i]->addUninitialized(1);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->addUninitialized(1);
}
/* ----------------------------------------------------------------------
@ -209,8 +176,8 @@
template<typename T>
void AssociativePointerArray<T>::addZeroElement()
{
for(int i=0;i<numElem_;i++)
content_[i]->addZero();
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->addZero();
}
/* ----------------------------------------------------------------------
@ -220,8 +187,8 @@
template<typename T>
void AssociativePointerArray<T>::deleteElement(int n)
{
for(int i=0;i<numElem_;i++)
content_[i]->del(n);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->del(n);
}
/* ----------------------------------------------------------------------
@ -231,8 +198,8 @@
template<typename T>
void AssociativePointerArray<T>::deleteForwardElement(int n,bool scale,bool translate,bool rotate)
{
for(int i=0;i<numElem_;i++)
content_[i]->delForward(n,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->delForward(n,scale,translate,rotate);
}
/* ----------------------------------------------------------------------
@ -242,8 +209,8 @@
template<typename T>
void AssociativePointerArray<T>::deleteRestartElement(int n,bool scale,bool translate,bool rotate)
{
for(int i=0;i<numElem_;i++)
content_[i]->delRestart(n,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->delRestart(n,scale,translate,rotate);
}
/* ----------------------------------------------------------------------
@ -253,8 +220,8 @@
template<typename T>
void AssociativePointerArray<T>::deleteRestartGlobal(bool scale,bool translate,bool rotate)
{
for(int i=0;i<numElem_;i++)
content_[i]->delRestart(scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->delRestart(scale,translate,rotate);
}
/* ----------------------------------------------------------------------
@ -264,27 +231,20 @@
template<typename T>
void AssociativePointerArray<T>::clearReverse(bool scale,bool translate,bool rotate)
{
for(int i=0;i<numElem_;i++)
content_[i]->clearReverse(scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->clearReverse(scale,translate,rotate);
}
/* ----------------------------------------------------------------------
id 2 index
has id
------------------------------------------------------------------------- */
template<typename T>
int AssociativePointerArray<T>::idToIndex(const char *_id)
bool AssociativePointerArray<T>::hasId(const char *_id)
{
for(int i=0;i<numElem_;i++)
if(content_[i]->matches_id(_id))
return i;
return -1;
}
content_iterator it = content_.find(_id);
return (it != content_.end());
template<typename T>
void AssociativePointerArray<T>::indexToId(int index, char *_id)
{
content_[index]->id(_id);
}
/* ----------------------------------------------------------------------
@ -294,17 +254,16 @@
template<typename T>
void AssociativePointerArray<T>::storeOrig(AssociativePointerArray &orig)
{
for(int i = 0; i < numElem_; i++)
orig.content_[i]->setFromContainer(content_[i]);
for(content_iterator it = orig.content_.begin(); it != orig.content_.end(); ++it)
it->second->setFromContainer(content_[it->first]);
}
template<typename T>
void AssociativePointerArray<T>::storeOrig(const char *_id, AssociativePointerArray &orig)
{
/*NL*/ //printf("storeOrig called for %s \n",_id);
for(int i = 0; i < numElem_; i++)
if(content_[i]->matches_id(_id))
orig.content_[i]->setFromContainer(content_[i]);
for(content_iterator it = orig.content_.begin(); it != orig.content_.end(); ++it)
if(content_[it->first]->matches_id(_id))
it->second->setFromContainer(content_[it->first]);
}
/* ----------------------------------------------------------------------
@ -314,24 +273,28 @@
template<typename T>
bool AssociativePointerArray<T>::reset(AssociativePointerArray &orig)
{
/*NL*/ //printf("numElem_ %d\n",numElem_);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->setFromContainer(orig.content_[it->first]);
for(int i = 0; i < numElem_; i++)
content_[i]->setFromContainer(orig.content_[i]);
return true;
return true;
}
template<typename T>
bool AssociativePointerArray<T>::reset(const char *_id, AssociativePointerArray &orig)
{
/*NL*/ //printf("numElem_ %d\n",numElem_);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
if(it->second->matches_id(_id))
it->second->setFromContainer(orig.content_[it->first]);
for(int i = 0; i < numElem_; i++)
if(content_[i]->matches_id(_id))
content_[i]->setFromContainer(orig.content_[i]);
return true;
}
return true;
template<typename T>
void AssociativePointerArray<T>::setToDefault(int n)
{
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
if(it->second->useDefault())
it->second->setToDefault(n);
}
/* ----------------------------------------------------------------------
@ -339,31 +302,31 @@
------------------------------------------------------------------------- */
template<typename T>
void AssociativePointerArray<T>::rotate(double *dQ)
void AssociativePointerArray<T>::rotate(const double *dQ)
{
for(int i = 0; i < numElem_; i++)
content_[i]->rotate(dQ);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->rotate(dQ);
}
template<typename T>
void AssociativePointerArray<T>::scale(double factor)
{
for(int i = 0; i < numElem_;i++)
content_[i]->scale(factor);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->scale(factor);
}
template<typename T>
void AssociativePointerArray<T>::move(double *delta)
void AssociativePointerArray<T>::move(const double *delta)
{
for(int i = 0; i < numElem_;i++)
content_[i]->move(delta);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->move(delta);
}
template<typename T>
void AssociativePointerArray<T>::moveElement(int n,double *delta)
void AssociativePointerArray<T>::moveElement(int n, const double *delta)
{
for(int i = 0; i < numElem_;i++)
content_[i]->moveElement(n,delta);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
it->second->moveElement(n,delta);
}
/* ----------------------------------------------------------------------
@ -374,8 +337,8 @@
int AssociativePointerArray<T>::bufSize(int operation,bool scale,bool translate,bool rotate) const
{
int buf_size = 0;
for(int i=0;i<numElem_;i++)
buf_size += getBasePointerByIndex(i)->bufSize(operation,scale,translate,rotate);
for(content_const_iterator it = content_.begin(); it != content_.end(); ++it)
buf_size += it->second->bufSize(operation,scale,translate,rotate);
return buf_size;
}
@ -383,8 +346,8 @@
int AssociativePointerArray<T>::pushToBuffer(double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nsend = 0;
for(int i=0;i<numElem_;i++)
nsend += getBasePointerByIndex(i)->pushToBuffer(&(buf[nsend]),operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nsend += it->second->pushToBuffer(&(buf[nsend]),operation,scale,translate,rotate);
return nsend;
}
@ -392,8 +355,8 @@
int AssociativePointerArray<T>::popFromBuffer(double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nrecv = 0;
for(int i=0;i<numElem_;i++)
nrecv += getBasePointerByIndex(i)->popFromBuffer(&(buf[nrecv]),operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nrecv += it->second->popFromBuffer(&(buf[nrecv]),operation,scale,translate,rotate);
return nrecv;
}
@ -405,8 +368,8 @@
int AssociativePointerArray<T>::elemListBufSize(int n,int operation,bool scale,bool translate,bool rotate)
{
int buf_size = 0;
for(int i=0;i<numElem_;i++)
buf_size += getBasePointerByIndex(i)->elemListBufSize(n,operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
buf_size += it->second->elemListBufSize(n,operation,scale,translate,rotate);
return buf_size;
}
@ -414,8 +377,8 @@
int AssociativePointerArray<T>::pushElemListToBuffer(int n, int *list, double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nsend = 0;
for(int i=0;i<numElem_;i++)
nsend += getBasePointerByIndex(i)->pushElemListToBuffer(n,list,&buf[nsend],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nsend += it->second->pushElemListToBuffer(n,list,&buf[nsend],operation,scale,translate,rotate);
return nsend;
}
@ -423,8 +386,8 @@
int AssociativePointerArray<T>::popElemListFromBuffer(int first, int n, double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nrecv = 0;
for(int i=0;i<numElem_;i++)
nrecv += getBasePointerByIndex(i)->popElemListFromBuffer(first,n,&buf[nrecv],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nrecv += it->second->popElemListFromBuffer(first,n,&buf[nrecv],operation,scale,translate,rotate);
return nrecv;
}
@ -432,8 +395,8 @@
int AssociativePointerArray<T>::pushElemListToBufferReverse(int first, int n, double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nrecv = 0;
for(int i=0;i<numElem_;i++)
nrecv += getBasePointerByIndex(i)->pushElemListToBufferReverse(first,n,&buf[nrecv],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nrecv += it->second->pushElemListToBufferReverse(first,n,&buf[nrecv],operation,scale,translate,rotate);
return nrecv;
}
@ -441,8 +404,8 @@
int AssociativePointerArray<T>::popElemListFromBufferReverse(int n, int *list, double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nsend = 0;
for(int i=0;i<numElem_;i++)
nsend += getBasePointerByIndex(i)->popElemListFromBufferReverse(n,list,&buf[nsend],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nsend += it->second->popElemListFromBufferReverse(n,list,&buf[nsend],operation,scale,translate,rotate);
return nsend;
}
@ -454,8 +417,8 @@
int AssociativePointerArray<T>::elemBufSize(int operation,bool scale,bool translate,bool rotate)
{
int buf_size = 0;
for(int i=0;i<numElem_;i++)
buf_size += getBasePointerByIndex(i)->elemBufSize(operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
buf_size += it->second->elemBufSize(operation,scale,translate,rotate);
return buf_size;
}
@ -463,8 +426,8 @@
int AssociativePointerArray<T>::pushElemToBuffer(int n, double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nsend = 0;
for(int i=0;i<numElem_;i++)
nsend += getBasePointerByIndex(i)->pushElemToBuffer(n,&buf[nsend],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nsend += it->second->pushElemToBuffer(n,&buf[nsend],operation,scale,translate,rotate);
return nsend;
}
@ -472,8 +435,8 @@
int AssociativePointerArray<T>::popElemFromBuffer(double *buf, int operation,bool scale,bool translate, bool rotate)
{
int nrecv = 0;
for(int i=0;i<numElem_;i++)
nrecv += getBasePointerByIndex(i)->popElemFromBuffer(&buf[nrecv],operation,scale,translate,rotate);
for(content_iterator it = content_.begin(); it != content_.end(); ++it)
nrecv += it->second->popElemFromBuffer(&buf[nrecv],operation,scale,translate,rotate);
return nrecv;
}

View File

@ -14,6 +14,7 @@
#include <string.h>
#include "bond.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "suffix.h"
@ -76,6 +77,8 @@ void Bond::n_granhistory(int nhist)
{
ngranhistory = nhist;
atom->n_bondhist = ngranhistory;
if (atom->nmax)
atom->avec->grow(atom->nmax); // make sure atom->bond_hist is created
//NP error check necessary because n_bondhist influences avec grow
//NP do not check this here since is already set before in AtomVecBondGran::settings()

View File

@ -35,6 +35,9 @@ namespace LAMMPS_NS
BoundingBox::BoundingBox()
: xLo(0.), xHi(0.), yLo(0.), yHi(0.), zLo(0.), zHi(0.), initGiven(false), dirty(true)
{}
BoundingBox::BoundingBox(double bounds[6])
: xLo(bounds[0]), xHi(bounds[1]), yLo(bounds[2]), yHi(bounds[3]), zLo(bounds[4]), zHi(bounds[5]), initGiven(true), dirty(true)
{}
BoundingBox::BoundingBox(double xLo_, double xHi_, double yLo_, double yHi_, double zLo_, double zHi_)
: xLo(xLo_), xHi(xHi_), yLo(yLo_), yHi(yHi_), zLo(zLo_), zHi(zHi_), initGiven(true), dirty(true)
{}

View File

@ -41,6 +41,7 @@ class BoundingBox
public:
BoundingBox();
explicit BoundingBox(double bounds[6]);
BoundingBox(double xLo, double xHi, double yLo, double yHi, double zLo, double zHi);
virtual ~BoundingBox();
@ -164,7 +165,7 @@ class BoundingBox
bool isInitialized()
{ return initGiven; }
bool isInside(double *p)
bool isInside(const double *p)
{
// check bbox
// test for >= and < as in Domain class

View File

@ -126,7 +126,7 @@ void CfdDatacoupling::grow_()
pull property from other code
------------------------------------------------------------------------- */
void CfdDatacoupling::pull(const char *name, const char *type, void *&, const char *)
void CfdDatacoupling::pull(const char *name, const char *type, void *&, const char *, int iworld)
{
// for MPI this is called by the library interface
// check if the requested property was registered by a LIGGGHTS model
@ -164,7 +164,7 @@ void CfdDatacoupling::pull(const char *name, const char *type, void *&, const ch
push property to other code
------------------------------------------------------------------------- */
void CfdDatacoupling::push(const char *name, const char *type, void *&, const char *)
void CfdDatacoupling::push(const char *name, const char *type, void *&, const char *, int iworld)
{
// for MPI this is called by the library interface
// check if the requested property was registered by a LIGGGHTS model

View File

@ -37,8 +37,8 @@ class CfdDatacoupling : protected Pointers {
void add_pull_property(const char*, const char*);
void add_push_property(const char*, const char*);
virtual void pull(const char *name, const char *type, void *&ptr, const char *datatype);
virtual void push(const char *name, const char *type, void *&ptr, const char *datatype);
virtual void pull(const char *name, const char *type, void *&ptr, const char *datatype, int iworld = 0);
virtual void push(const char *name, const char *type, void *&ptr, const char *datatype, int iworld = 0);
virtual void allocate_external(int **&data, int len2,int len1,int initvalue);
virtual void allocate_external(double **&data, int len2,int len1,double initvalue);

View File

@ -113,7 +113,7 @@ void CfdDatacouplingFile::exchange()
//NP called with from = NULL
void CfdDatacouplingFile::pull(const char *name, const char *type, void *&from, const char *datatype)
void CfdDatacouplingFile::pull(const char *name, const char *type, void *&from, const char *datatype, int iworld)
{
CfdDatacoupling::pull(name,type,from,datatype);
@ -152,7 +152,7 @@ void CfdDatacouplingFile::pull(const char *name, const char *type, void *&from,
//NP called with to = NULL
void CfdDatacouplingFile::push(const char *name, const char *type, void *&to, const char *datatype)
void CfdDatacouplingFile::push(const char *name, const char *type, void *&to, const char *datatype, int iworld)
{
CfdDatacoupling::push(name,type,to,datatype);

View File

@ -38,8 +38,8 @@ class CfdDatacouplingFile : public CfdDatacoupling {
~CfdDatacouplingFile();
friend class FixTempFromFile;
void pull(const char *, const char *, void *&, const char *);
void push(const char *, const char *, void *&, const char *);
void pull(const char *, const char *, void *&, const char *, int iworld = 0);
void push(const char *, const char *, void *&, const char *, int iworld = 0);
virtual void post_create();
void exchange();

View File

@ -74,27 +74,27 @@ void CfdDatacouplingMPI::exchange()
/* ---------------------------------------------------------------------- */
void CfdDatacouplingMPI::pull(const char *name,const char *type,void *&from,const char *datatype)
void CfdDatacouplingMPI::pull(const char *name,const char *type,void *&from,const char *datatype,int iworld)
{
CfdDatacoupling::pull(name,type,from,datatype);
if(strcmp(datatype,"double") == 0)
pull_mpi<double>(name,type,from);
pull_mpi<double>(name,type,from,iworld);
else if(strcmp(datatype,"int") == 0)
pull_mpi<int>(name,type,from);
pull_mpi<int>(name,type,from,iworld);
else error->one(FLERR,"Illegal call to CfdDatacouplingMPI::pull, valid datatypes are 'int' and double'");
}
/* ---------------------------------------------------------------------- */
void CfdDatacouplingMPI::push(const char *name,const char *type,void *&to,const char *datatype)
void CfdDatacouplingMPI::push(const char *name,const char *type,void *&to,const char *datatype,int iworld)
{
CfdDatacoupling::push(name,type,to,datatype);
if(strcmp(datatype,"double") == 0)
push_mpi<double>(name,type,to);
push_mpi<double>(name,type,to,iworld);
else if(strcmp(datatype,"int") == 0)
push_mpi<int>(name,type,to);
push_mpi<int>(name,type,to,iworld);
else error->one(FLERR,"Illegal call to CfdDatacouplingMPI::pull, valid datatypes are 'int' and double'");
}

View File

@ -28,10 +28,13 @@
#ifndef LMP_CFD_DATACOUPLING_MPI_H
#define LMP_CFD_DATACOUPLING_MPI_H
#define MULTI_PARTITION_CFD
#include "cfd_datacoupling.h"
#include "multisphere_parallel.h"
#include "error.h"
#include "properties.h"
#include "universe.h"
#include <mpi.h>
namespace LAMMPS_NS {
@ -43,11 +46,11 @@ class CfdDatacouplingMPI : public CfdDatacoupling {
void exchange();
virtual void pull(const char *name, const char *type, void *&ptr, const char *datatype);
virtual void push(const char *name, const char *type, void *&ptr, const char *datatype);
virtual void pull(const char *name, const char *type, void *&ptr, const char *datatype,int iworld=0);
virtual void push(const char *name, const char *type, void *&ptr, const char *datatype,int iworld=0);
template <typename T> void pull_mpi(const char *,const char *,void *&);
template <typename T> void push_mpi(const char *,const char *,void *&);
template <typename T> void pull_mpi(const char *,const char *,void *&,int iworld=0);
template <typename T> void push_mpi(const char *,const char *,void *&,int iworld=0);
virtual bool error_push()
{ return false;}
@ -74,29 +77,57 @@ class CfdDatacouplingMPI : public CfdDatacoupling {
//NP ptr to is local, ptr from is global
template <typename T>
void CfdDatacouplingMPI::pull_mpi(const char *name,const char *type,void *&from)
void CfdDatacouplingMPI::pull_mpi(const char *name,const char *type,void *&from,int iworld)
{
int len1 = -1, len2 = -1, m;
// len1 = atom->tag_max(); except for scalar-global, vector-global, matrix-global
// get reference where to write the data
void * to = find_pull_property(name,type,len1,len2);
#ifdef MULTI_PARTITION_CFD
if (universe->existflag == 0 || iworld == universe->iworld)
#endif
{
if (atom->nlocal && (!to || len1 < 0 || len2 < 0))
{
if(screen) fprintf(screen,"LIGGGHTS could not find property %s to write data from calling program to.\n",name);
lmp->error->one(FLERR,"This is fatal");
}
}
int total_len = len1*len2;
#ifdef MULTI_PARTITION_CFD
if (universe->existflag == 1) // enforce size of requested world
MPI_Bcast(&total_len, 1, MPI_INT, universe->root_proc[iworld], universe->uworld);
#endif
// return if no data to transmit
if(len1*len2 < 1) return;
if(total_len < 1) return;
// check memory allocation
T* allred = check_grow<T>(len1*len2);
T* allred = check_grow<T>(total_len);
// perform allreduce on incoming data
T **from_t = (T**)from;
MPI_Allreduce(&(from_t[0][0]),&(allred[0]),len1*len2,mpi_type_dc<T>(),MPI_SUM,world);
#ifndef MULTI_PARTITION_CFD
MPI_Allreduce(&(from_t[0][0]), &(allred[0]), total_len, mpi_type_dc<T>(), MPI_SUM, world);
#else
//MPI_Reduce + MPI_Bcast
// int MPI_Reduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm)
// int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm)
//reduce from all procs in universe onto root proc of iworld
MPI_Reduce (&(from_t[0][0]), &(allred[0]), total_len, mpi_type_dc<T>(), MPI_SUM, universe->root_proc[iworld], universe->uworld);
//bcast to iworld
if (iworld == universe->iworld) {
MPI_Bcast (&(allred[0]), total_len, mpi_type_dc<T>(), 0, world);
}
#endif
#ifdef MULTI_PARTITION_CFD
if(iworld == universe->iworld) // only copy to requested world
#endif
{
// copy data - loops over max # global atoms, bodies
if(strcmp(type,"scalar-atom") == 0)
{
@ -142,6 +173,7 @@ void CfdDatacouplingMPI::pull_mpi(const char *name,const char *type,void *&from)
to_t[i][j] = allred[i*len2 + j];
}
else error->one(FLERR,"Illegal data type in CfdDatacouplingMPI::pull");
}
}
/* ---------------------------------------------------------------------- */
@ -150,7 +182,7 @@ void CfdDatacouplingMPI::pull_mpi(const char *name,const char *type,void *&from)
//NP len1 is global # of datums (max tag)
template <typename T>
void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to)
void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to,int iworld)
{
int len1 = -1, len2 = -1, id;
@ -164,6 +196,9 @@ void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to)
// get reference where to write the data
void * from = find_push_property(name,type,len1,len2);
#ifdef MULTI_PARTITION_CFD
if (universe->existflag == 0 || iworld == universe->iworld)
#endif
if (atom->nlocal && (!from || len1 < 0 || len2 < 0))
{
/*NL*/ //if(screen) fprintf(screen,"nlocal %d, len1 %d lens2 %d\n",atom->nlocal,len1,len2);
@ -171,15 +206,24 @@ void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to)
lmp->error->one(FLERR,"This is fatal");
}
int total_len = len1*len2;
#ifdef MULTI_PARTITION_CFD
if (universe->existflag == 1) // enforce size of requested world
MPI_Bcast(&total_len, 1, MPI_INT, universe->root_proc[iworld], universe->uworld);
#endif
// return if no data to transmit
if(len1*len2 < 1) return;
if(total_len < 1) return;
// check memory allocation
T * allred = check_grow<T>(len1*len2);
T * allred = check_grow<T>(total_len);
// zeroize before using allreduce
vectorZeroizeN(allred,len1*len2);
vectorZeroizeN(allred,total_len);
#ifdef MULTI_PARTITION_CFD
if(iworld == universe->iworld) // only copy from requested world
#endif
{
// copy data - loop local # atoms, bodies
if(strcmp(type,"scalar-atom") == 0)
{
@ -225,7 +269,7 @@ void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to)
/*NL*/ //if(screen) fprintf(screen,"id %d from %d %d: %f\n",id ,i,j,from_t[i][j]);
}
}
/*NL*/ //if (screen) printVecN(screen,"allred",allred,len1*len2);
/*NL*/ //if (screen) printVecN(screen,"allred",allred,total_len);
}
else if(strcmp(type,"scalar-global") == 0 || strcmp(type,"vector-global") == 0 || strcmp(type,"matrix-global") == 0)
{
@ -235,10 +279,24 @@ void CfdDatacouplingMPI::push_mpi(const char *name,const char *type,void *&to)
allred[i*len2 + j] = from_t[i][j];
}
else error->one(FLERR,"Illegal data type in CfdDatacouplingMPI::pull");
}
// perform allreduce on outgoing data
T **to_t = (T**)to;
MPI_Allreduce(&(allred[0]),&(to_t[0][0]),len1*len2,mpi_type_dc<T>(),MPI_SUM,world);
#ifndef MULTI_PARTITION_CFD
MPI_Allreduce(&(allred[0]),&(to_t[0][0]),total_len,mpi_type_dc<T>(),MPI_SUM,world);
#else
//MPI_Reduce + MPI_Bcast
// int MPI_Reduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm)
// int MPI_Bcast( void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm)
// reduce from request world onto root of requested world
if (iworld == universe->iworld) {
MPI_Reduce(&(allred[0]), &(to_t[0][0]), total_len, mpi_type_dc<T>(), MPI_SUM, 0, world);
}
// bcast from root of requested world to all procs in universe
MPI_Bcast(&(to_t[0][0]), total_len, mpi_type_dc<T>(), universe->root_proc[iworld], universe->uworld);
#endif
}
/* ---------------------------------------------------------------------- */

View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <string.h>
#include "compute_com_molecule.h"
#include "atom.h"
#include "update.h"

View File

@ -0,0 +1,60 @@
/* ----------------------------------------------------------------------
LIGGGHTS - LAMMPS Improved for General Granular and Granular Heat
Transfer Simulations
Copyright 2015- JKU Linz
LIGGGHTS is based on LAMMPS
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors:
Daniel Queteschiner (JKU Linz)
------------------------------------------------------------------------- */
#ifndef LMP_CONST_PARTICLE_TEMPLATE_SPHERE_H
#define LMP_CONST_PARTICLE_TEMPLATE_SPHERE_H
#include <cmath>
#include <map>
#include "math_const.h"
namespace LAMMPS_NS {
class ConstantParticleTemplateSphere
{
public:
ConstantParticleTemplateSphere() :
radius_(0.), density_(0.), mass_(0), atomtype_(0) {}
ConstantParticleTemplateSphere(double radius, double mass, int atomtype) :
radius_(radius), mass_(mass), atomtype_(atomtype) { density_ = mass/((4./3.)*MathConst::MY_PI*radius*radius*radius); }
double radius_, density_, mass_;
int atomtype_;
bool operator<(const ConstantParticleTemplateSphere &other) const
{
//return radius < other.radius; // sort smallest to largest radius
if(radius_ > other.radius_) return true; // sort largest to smallest radius
else if(radius_ < other.radius_) return false;
if(mass_ > other.mass_) return true; // sort largest to smallest mass
else if(mass_ < other.mass_) return false;
if(atomtype_ > other.atomtype_) return true;
/*else if(atomtype_ < other.atomtype_) return false;*/
return false; // always return false for equivalent templates
}
};
typedef std::map<ConstantParticleTemplateSphere,double> DiscreteParticleDistribution; // template, #particles
}
#endif

View File

@ -78,9 +78,9 @@ namespace LAMMPS_NS
virtual bool setFromContainer(ContainerBase *cont) = 0;
virtual void scale(double factor) = 0;
virtual void move(double *dx) = 0;
virtual void moveElement(int i,double *dx) = 0;
virtual void rotate(double *dQ) = 0;
virtual void move(const double *dx) = 0;
virtual void moveElement(int i, const double *dx) = 0;
virtual void rotate(const double *dQ) = 0;
virtual void setToDefault(int n) = 0;

View File

@ -117,7 +117,7 @@ using namespace LAMMPS_NS;
rotate all properties, applies to vector and multivector only
------------------------------------------------------------------------- */
void CustomValueTracker::rotate(double *totalQ,double *dQ)
void CustomValueTracker::rotate(const double *totalQ, const double *dQ)
{
/*NL*/ //if (screen) printVec4D(screen,"totalQ",totalQ);
/*NL*/ //if (screen) printVec4D(screen,"dQ",dQ);
@ -127,7 +127,7 @@ using namespace LAMMPS_NS;
globalProperties_.rotate(totalQ);
}
void CustomValueTracker::rotate(double *dQ)
void CustomValueTracker::rotate(const double *dQ)
{
//NP this handles owned and ghost elements
elementProperties_.rotate(dQ);
@ -149,14 +149,14 @@ using namespace LAMMPS_NS;
move all properties
------------------------------------------------------------------------- */
void CustomValueTracker::move(double *vecTotal, double *vecIncremental)
void CustomValueTracker::move(const double *vecTotal, const double *vecIncremental)
{
//NP this handles owned and ghost elements
elementProperties_.move(vecIncremental);
globalProperties_.move(vecTotal);
}
void CustomValueTracker::move(double *vecIncremental)
void CustomValueTracker::move(const double *vecIncremental)
{
//NP this handles owned and ghost elements
elementProperties_.move(vecIncremental);

View File

@ -50,9 +50,10 @@ namespace LAMMPS_NS
T* getElementProperty(const char *_id);
inline ContainerBase* getElementPropertyBase(const char *_id);
inline ContainerBase* getElementPropertyBase(int i);
inline int getElementPropertyIndex(const char *_id);
inline bool hasElementProperty(const char *_id);
inline void setElementPropertyToDefault(int n);
template<typename T, typename U>
void setElementProperty(const char *_id, U def);
@ -92,11 +93,11 @@ namespace LAMMPS_NS
inline void storeGlobalPropOrig(const char *_id);
inline void resetGlobalPropToOrig(const char *_id);
inline void moveElement(int i, double *delta);
void move(double *vecTotal, double *vecIncremental);
void move(double *vecIncremental);
void rotate(double *totalQ, double *dQ);
void rotate(double *dQ);
inline void moveElement(int i, const double *delta);
void move(const double *vecTotal, const double *vecIncremental);
void move(const double *vecIncremental);
void rotate(const double *totalQ, const double *dQ);
void rotate(const double *dQ);
void scale(double factor);
// buffer operations

View File

@ -131,16 +131,17 @@
return elementProperties_.getBasePointerById(_id);
}
inline ContainerBase* CustomValueTracker::getElementPropertyBase(int i)
inline bool CustomValueTracker::hasElementProperty(const char *_id)
{
return elementProperties_.getBasePointerByIndex(i);
return elementProperties_.hasId(_id);
}
inline int CustomValueTracker::getElementPropertyIndex(const char *_id)
inline void CustomValueTracker::setElementPropertyToDefault(int n)
{
return elementProperties_.idToIndex(_id);
elementProperties_.setToDefault(n);
}
template<typename T>
T* CustomValueTracker::getGlobalProperty(const char *_id)
{
@ -259,7 +260,7 @@
move element i
------------------------------------------------------------------------- */
void CustomValueTracker::moveElement(int i, double *delta)
void CustomValueTracker::moveElement(int i, const double *delta)
{
elementProperties_.moveElement(i,delta);
}

View File

@ -166,19 +166,19 @@ class Domain : protected Pointers {
}
//NP modified C.K in domain_I.h
int is_in_domain(double* pos); //NP modified C.K.
int is_in_subdomain(double* pos); //NP modified C.K.
int is_in_extended_subdomain(double* pos); //NP modified C.K.
double dist_subbox_borders(double* pos); //NP modified C.K.
int is_in_domain(const double* pos); //NP modified C.K.
int is_in_subdomain(const double* pos); //NP modified C.K.
int is_in_extended_subdomain(const double* pos); //NP modified C.K.
double dist_subbox_borders(const double* pos); //NP modified C.K.
void min_subbox_extent(double &min_extent,int &dim); //NP modified C.K.
int is_periodic_ghost(int i); //NP modified C.K.
bool is_owned_or_first_ghost(int i); //NP modified C.K.
//NP modified C.K for wedge
virtual int is_in_domain_wedge(double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual int is_in_subdomain_wedge(double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual int is_in_extended_subdomain_wedge(double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual double dist_subbox_borders_wedge(double* pos) { UNUSED(pos); return 0.; } //NP modified C.K.
virtual int is_in_domain_wedge(const double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual int is_in_subdomain_wedge(const double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual int is_in_extended_subdomain_wedge(const double* pos) { UNUSED(pos); return 0; } //NP modified C.K.
virtual double dist_subbox_borders_wedge(const double* pos) { UNUSED(pos); return 0.; } //NP modified C.K.
virtual int is_periodic_ghost_wedge(int i) { UNUSED(i); return 0;} //NP modified C.K.
bool is_wedge; //NP modified C.K.

View File

@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
inlined for performance
------------------------------------------------------------------------- */
inline int Domain::is_in_domain(double* pos) //NP modified C.K.
inline int Domain::is_in_domain(const double* pos) //NP modified C.K.
{
if(is_wedge)
return is_in_domain_wedge(pos);
@ -45,7 +45,7 @@ inline int Domain::is_in_domain(double* pos) //NP modified C.K.
return 0;
}
inline int Domain::is_in_subdomain(double* pos) //NP modified C.K.
inline int Domain::is_in_subdomain(const double* pos) //NP modified C.K.
{
if(is_wedge)
return is_in_subdomain_wedge(pos);
@ -77,7 +77,7 @@ inline int Domain::is_in_subdomain(double* pos) //NP modified C.K.
return 0;
}
inline int Domain::is_in_extended_subdomain(double* pos) //NP modified C.K.
inline int Domain::is_in_extended_subdomain(const double* pos) //NP modified C.K.
{
if(is_wedge)
return is_in_extended_subdomain_wedge(pos);
@ -114,7 +114,7 @@ inline int Domain::is_in_extended_subdomain(double* pos) //NP modified C.K.
check distance from borders of subbox
------------------------------------------------------------------------- */
inline double Domain::dist_subbox_borders(double* pos) //NP modified C.K.
inline double Domain::dist_subbox_borders(const double* pos) //NP modified C.K.
{
if(is_wedge)
return dist_subbox_borders_wedge(pos);

View File

@ -358,8 +358,8 @@ void DumpATOMVTK::setFileCurrent() {
else {
char bif[8],pad[16];
strcpy(bif,BIGINT_FORMAT);
sprintf(pad,"%%s%%0%d_%%d%s%%s",padflag,&bif[1]);
sprintf(filecurrent,pad,filename,comm->me,update->ntimestep,ptr+1);
sprintf(pad,"%%s%%0%d%s%%s",padflag,&bif[1]);
sprintf(filecurrent,pad,filename,update->ntimestep,ptr+1);
}
*ptr = '*';
}

View File

@ -47,6 +47,37 @@ DumpEulerVTK::DumpEulerVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg,
if (narg < 5)
error->all(FLERR,"Illegal dump euler/vtk command");
fix_euler_name_ = NULL;
cell_center=true;
int iarg = 5;
while (iarg < narg) {
if (strcmp(arg[iarg], "ave_euler") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR,"missing argument");
++iarg;
int n = strlen(arg[iarg]) + 1;
fix_euler_name_ = new char[n];
strcpy(fix_euler_name_,arg[iarg]);
++iarg;
} else if (strcmp(arg[iarg], "cell_center") == 0) {
if (iarg + 1 >= narg)
error->all(FLERR,"missing argument");
++iarg;
if (strcmp(arg[iarg], "yes") == 0) {
cell_center = true;
} else if (strcmp(arg[iarg], "no") == 0) {
cell_center = false;
} else {
error->all(FLERR,"expecting 'yes' or 'no' after cell_center");;
}
++iarg;
} else {
error->all(FLERR,"illegal argument");;
}
}
// CURRENTLY ONLY PROC 0 writes
format_default = NULL;
@ -56,13 +87,20 @@ DumpEulerVTK::DumpEulerVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg,
DumpEulerVTK::~DumpEulerVTK()
{
delete fix_euler_name_;
}
/* ---------------------------------------------------------------------- */
void DumpEulerVTK::init_style()
{
fix_euler_ = static_cast<FixAveEuler*>(modify->find_fix_style("ave/euler",0));
if (fix_euler_name_) {
fix_euler_ = static_cast<FixAveEuler*>(modify->find_fix_id(fix_euler_name_));
}
if(!fix_euler_) {
if (fix_euler_name_) error->warning(FLERR, "dump euler/vtk failed to find named fix ave/euler*");
fix_euler_ = static_cast<FixAveEuler*>(modify->find_fix_style("ave/euler",0));
}
if(!fix_euler_)
error->all(FLERR,"Illegal dump euler/vtk command, need a fix ave/euler");
@ -74,12 +112,27 @@ void DumpEulerVTK::init_style()
error->all(FLERR,"Your 'dump euler/vtk' command is writing one file per processor, where all the files contain the same data");
// if (domain->triclinic == 1)
// error->all(FLERR,"Can not dump VTK files for triclinic box");
// error->all(FLERR,"Cannot dump VTK files for triclinic box");
if (binary)
error->all(FLERR,"Can not dump VTK files in binary mode");
error->all(FLERR,"Cannot dump VTK files in binary mode");
// node center (3), av vel (3), volume fraction, radius, stress
size_one = 9;
if (cell_center) {
// node center (3), av vel (3), volume fraction, radius, pressure, normal and shear stresses (6)
size_one = 15;
} else {
if (fix_euler_->is_parallel())
error->all(FLERR,"Cannot dump cell information for fix ave/euler running parallel");
if (strcmp(fix_euler_->style,"ave/euler") == 0) {
// boxlo (3), av vel (3), volume fraction, radius, pressure, normal and shear stresses (6)
size_one = 15;
} else if (strcmp(fix_euler_->style,"ave/euler/region") == 0) {
// node points (24), av vel (3), volume fraction, radius, pressure, normal and shear stresses (6)
size_one = 36;
} else {
error->all(FLERR,"Unable to dump cells to VTK files");
}
}
delete [] format;
}
@ -127,9 +180,25 @@ void DumpEulerVTK::pack(int *ids)
for(int i = 0; i < ncells; i++)
{
buf[m++] = fix_euler_->cell_center(i,0);
buf[m++] = fix_euler_->cell_center(i,1);
buf[m++] = fix_euler_->cell_center(i,2);
if (cell_center) {
buf[m++] = fix_euler_->cell_center(i,0);
buf[m++] = fix_euler_->cell_center(i,1);
buf[m++] = fix_euler_->cell_center(i,2);
} else {
if (strcmp(fix_euler_->style,"ave/euler") == 0) {
buf[m++] = fix_euler_->cell_center(i,0)-0.5*fix_euler_->cell_size(0);
buf[m++] = fix_euler_->cell_center(i,1)-0.5*fix_euler_->cell_size(1);
buf[m++] = fix_euler_->cell_center(i,2)-0.5*fix_euler_->cell_size(2);
} else if (strcmp(fix_euler_->style,"ave/euler/region") == 0) {
double points[24];
fix_euler_->cell_points(i,points);
for (int j=0; j < 24; ++j) {
buf[m++] = points[j];
}
} else {
error->all(FLERR,"Unable to dump cells to VTK files");
}
}
buf[m++] = fix_euler_->cell_v_av(i,0);
buf[m++] = fix_euler_->cell_v_av(i,1);
@ -138,6 +207,13 @@ void DumpEulerVTK::pack(int *ids)
buf[m++] = fix_euler_->cell_vol_fr(i);
buf[m++] = fix_euler_->cell_radius(i);
buf[m++] = fix_euler_->cell_pressure(i);
buf[m++] = fix_euler_->cell_stress(i,0);
buf[m++] = fix_euler_->cell_stress(i,1);
buf[m++] = fix_euler_->cell_stress(i,2);
buf[m++] = fix_euler_->cell_stress(i,3);
buf[m++] = fix_euler_->cell_stress(i,4);
buf[m++] = fix_euler_->cell_stress(i,5);
}
return ;
}
@ -178,25 +254,106 @@ void DumpEulerVTK::write_data_ascii(int n, double *mybuf)
/*NL*///error->one(FLERR,"end");
// write point data
fprintf(fp,"DATASET POLYDATA\nPOINTS %d float\n",n);
m = 0;
buf_pos = 0;
for (int i = 0; i < n; i++)
{
if (cell_center) {
fprintf(fp,"DATASET POLYDATA\nPOINTS %d float\n",n);
for (int i = 0; i < n; i++)
{
fprintf(fp,"%f %f %f\n",mybuf[m],mybuf[m+1],mybuf[m+2]);
m += size_one ;
}
buf_pos += 3;
}
buf_pos += 3;
// write polygon data
fprintf(fp,"VERTICES %d %d\n",n,2*n);
for (int i = 0; i < n; i++)
{
fprintf(fp,"%d %d\n",1,i);
}
// write polygon data
fprintf(fp,"VERTICES %d %d\n",n,2*n);
for (int i = 0; i < n; i++)
{
fprintf(fp,"%d %d\n",1,i);
}
// write point data header
fprintf(fp,"POINT_DATA %d\n",n);
// write point data header
fprintf(fp,"POINT_DATA %d\n",n);
} else {
if (strcmp(fix_euler_->style,"ave/euler") == 0) {
// dump a rectilinear grid of cuboidal cell
int nx = fix_euler_->ncells(0);
int ny = fix_euler_->ncells(1);
int nz = fix_euler_->ncells(2);
if (nx*ny*nz != n)
error->all(FLERR,"number of cells does not match nx*ny*nz != n");
fprintf(fp,"DATASET RECTILINEAR_GRID\nDIMENSIONS %d %d %d\n",nx+1,ny+1,nz+1);
// x coordinates
m = 0;
fprintf(fp,"X_COORDINATES %d float\n", nx+1);
for (int i = 0; i < nx; i++) {
fprintf(fp,"%f ",mybuf[m]);
m += size_one ;
}
fprintf(fp,"%f\n", mybuf[m-size_one]+fix_euler_->cell_size(0));
// y coordinates
m = 1;
fprintf(fp,"Y_COORDINATES %d float\n", ny+1);
for (int i = 0; i < ny; i++) {
fprintf(fp,"%f ",mybuf[m]);
m += nx*size_one ;
}
fprintf(fp,"%f\n", mybuf[m-nx*size_one]+fix_euler_->cell_size(1));
// z coordinates
m = 2;
fprintf(fp,"Z_COORDINATES %d float\n", nz+1);
for (int i = 0; i < nz; i++) {
fprintf(fp,"%f ",mybuf[m]);
m += nx*ny*size_one ;
}
fprintf(fp,"%f\n", mybuf[m-nx*ny*size_one]+fix_euler_->cell_size(2));
buf_pos += 3;
// write cell data header
fprintf(fp,"CELL_DATA %d\n",n);
} else if (strcmp(fix_euler_->style,"ave/euler/region") == 0) {
// dump an unstructured grid of hexahedral cells
fprintf(fp,"DATASET UNSTRUCTURED_GRID\nPOINTS %d float\n",n*8);
for (int i = 0; i < n; i++) {
fprintf(fp,"%f %f %f %f %f %f\n%f %f %f %f %f %f\n"
"%f %f %f %f %f %f\n%f %f %f %f %f %f\n",
mybuf[m], mybuf[m+1], mybuf[m+2],
mybuf[m+3], mybuf[m+4], mybuf[m+5],
mybuf[m+6], mybuf[m+7], mybuf[m+8],
mybuf[m+9], mybuf[m+10],mybuf[m+11],
mybuf[m+12],mybuf[m+13],mybuf[m+14],
mybuf[m+15],mybuf[m+16],mybuf[m+17],
mybuf[m+18],mybuf[m+19],mybuf[m+20],
mybuf[m+21],mybuf[m+22],mybuf[m+23]);
m += size_one ;
}
fprintf(fp,"CELLS %d %d\n",n,n*9);
for (int i = 0; i < n; i++) {
fprintf(fp,"8 %d %d %d %d %d %d %d %d\n",
i*8,i*8+1,i*8+2,i*8+3,i*8+4,i*8+5,i*8+6,i*8+7);
}
fprintf(fp,"CELL_TYPES %d\n",n);
for (int i = 0; i < n; i++) {
fprintf(fp,"12\n"); // hexahedron = 12
}
buf_pos += 24;
// write cell data header
fprintf(fp,"CELL_DATA %d\n",n);
} else {
error->all(FLERR,"Unable to dump cells to VTK files");
}
}
// write cell data
@ -236,6 +393,27 @@ void DumpEulerVTK::write_data_ascii(int n, double *mybuf)
}
buf_pos++;
fprintf(fp,"VECTORS sigma float\n");
m = buf_pos;
for (int i = 0; i < n; i++)
{
// normal stresses
fprintf(fp,"%f %f %f\n",mybuf[m],mybuf[m+1],mybuf[m+2]);
m += size_one;
}
buf_pos += 3;
fprintf(fp,"VECTORS tau float\n");
m = buf_pos;
for (int i = 0; i < n; i++)
{
// shear stresses
// note: internal order of stresses is (xx, yy, zz, xy, xz, yz)
// Voigt notation is (xx, yy, zz, yz, xz, xy)
fprintf(fp,"%f %f %f\n",mybuf[m],mybuf[m+1],mybuf[m+2]);
m += size_one;
}
buf_pos += 3;
// footer not needed
// if would be needed, would do like in dump stl
}

View File

@ -58,6 +58,8 @@ class DumpEulerVTK : public Dump {
// buffer for data from all procs
int n_all_, n_all_max_;
double *buf_all_;
char *fix_euler_name_;
bool cell_center;
};

View File

@ -81,6 +81,9 @@ DumpMeshVTK::DumpMeshVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, ar
f_node_(0),
T_(0),
min_active_edge_dist_(0),
int_scalar_containers_(0),
int_scalar_container_names_(0),
n_int_scalar_containers_(0),
scalar_containers_(0),
scalar_container_names_(0),
n_scalar_containers_(0),
@ -249,19 +252,24 @@ DumpMeshVTK::DumpMeshVTK(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, ar
min_active_edge_dist_ = new ScalarContainer<double>*[nMesh_];
//NP allocate container_bases_
int_scalar_containers_ = new ScalarContainer<int>**[n_container_bases_];
int_scalar_container_names_ = new char*[n_container_bases_];
scalar_containers_ = new ScalarContainer<double>**[n_container_bases_];
scalar_container_names_ = new char*[n_container_bases_];
vector_containers_ = new VectorContainer<double,3>**[n_container_bases_];
vector_container_names_ = new char*[n_container_bases_];
for(int i = 0; i < n_container_bases_; i++)
{
int_scalar_containers_[i] = new ScalarContainer<int>*[nMesh_];
scalar_containers_[i] = new ScalarContainer<double>*[nMesh_];
vector_containers_[i] = new VectorContainer<double,3>*[nMesh_];
int_scalar_container_names_[i] = new char[200];
scalar_container_names_[i] = new char[200];
vector_container_names_[i] = new char[200];
for(int j = 0; j < nMesh_; j++)
{
int_scalar_containers_[i][j] = 0;
scalar_containers_[i][j] = 0;
vector_containers_[i][j] = 0;
}
@ -293,13 +301,17 @@ DumpMeshVTK::~DumpMeshVTK()
for(int i = 0; i < n_container_bases_; i++)
{
delete [] int_scalar_containers_[i];
delete [] scalar_containers_[i];
delete [] vector_containers_[i];
delete [] int_scalar_container_names_[i];
delete [] scalar_container_names_[i];
delete [] vector_container_names_[i];
}
delete [] int_scalar_containers_;
delete [] scalar_containers_;
delete [] vector_containers_;
delete [] int_scalar_container_names_;
delete [] scalar_container_names_;
delete [] vector_container_names_;
}
@ -361,6 +373,7 @@ void DumpMeshVTK::init_style()
//NP get refs so can compute size
getGeneralRefs();
size_one += n_int_scalar_containers_*1;
size_one += n_scalar_containers_*1;
size_one += n_vector_containers_*3;
@ -478,19 +491,28 @@ void DumpMeshVTK::getGeneralRefs()
//NP general implementation
//NP go through list of fields
n_int_scalar_containers_ = 0;
n_scalar_containers_ = 0;
n_vector_containers_ = 0;
char cid[200];
for(int ib = 0; ib < n_container_bases_; ib++)
{
bool found_scalar = false, found_vector = false;
bool found_int_scalar = false, found_scalar = false, found_vector = false;
/*NL*/ //if (screen) fprintf(screen,"looking for container %s\n",container_args_[ib]);
for(int i = 0; i < nMesh_; i++)
{
if(meshList_[i]->prop().getElementProperty<ScalarContainer<double> >(container_args_[ib]))
if(meshList_[i]->prop().getElementProperty<ScalarContainer<int> >(container_args_[ib]))
{
/*NL*/ //fprintf(screen,"mesh %s: prop %s found\n",meshList_[i]->mesh_id(),container_args_[ib]);
found_int_scalar = true;
int_scalar_containers_[n_int_scalar_containers_][i] = meshList_[i]->prop().getElementProperty<ScalarContainer<int> >(container_args_[ib]);
int_scalar_containers_[n_int_scalar_containers_][i]->id(cid);
strcpy(int_scalar_container_names_[n_int_scalar_containers_],cid);
}
else if(meshList_[i]->prop().getElementProperty<ScalarContainer<double> >(container_args_[ib]))
{
/*NL*/ //if (screen) fprintf(screen,"mesh %s: prop %s found\n",meshList_[i]->mesh_id(),container_args_[ib]);
found_scalar = true;
@ -498,9 +520,7 @@ void DumpMeshVTK::getGeneralRefs()
scalar_containers_[n_scalar_containers_][i]->id(cid);
strcpy(scalar_container_names_[n_scalar_containers_],cid);
}
/*NL*/ //else if (screen) fprintf(screen,"mesh %s: prop %s NOT found\n",meshList_[i]->mesh_id(),container_args_[ib]);
if(meshList_[i]->prop().getElementProperty<VectorContainer<double,3> >(container_args_[ib]))
else if(meshList_[i]->prop().getElementProperty<VectorContainer<double,3> >(container_args_[ib]))
{
found_vector = true;
vector_containers_[n_vector_containers_][i] = meshList_[i]->prop().getElementProperty<VectorContainer<double,3> >(container_args_[ib]);
@ -509,12 +529,14 @@ void DumpMeshVTK::getGeneralRefs()
}
}
if(found_int_scalar)
n_int_scalar_containers_++;
if(found_scalar)
n_scalar_containers_++;
if(found_vector)
n_vector_containers_++;
if(!found_scalar && !found_vector)
if(!found_int_scalar && !found_scalar && !found_vector)
error->all(FLERR,"Illegal dump mesh/vtk command, unknown keyword or mesh");
}
@ -528,6 +550,7 @@ void DumpMeshVTK::pack(int *ids)
int m = 0;
double node[3];
TriMesh *mesh;
ScalarContainer<int> *cb0;
ScalarContainer<double> *cb1;
VectorContainer<double,3> *cb2;
@ -649,6 +672,12 @@ void DumpMeshVTK::pack(int *ids)
buf[m++] = min_active_edge_dist_[iMesh] ? min_active_edge_dist_[iMesh]->get(iTri) : 0.;
}
for(int ib = 0; ib < n_int_scalar_containers_; ib++)
{
cb0 = int_scalar_containers_[ib][iMesh];
buf[m++] = cb0 ? static_cast<double>((*cb0)(iTri)) : 0.;
}
for(int ib = 0; ib < n_scalar_containers_; ib++)
{
cb1 = scalar_containers_[ib][iMesh];
@ -1169,6 +1198,21 @@ void DumpMeshVTK::write_data_ascii_face(int n, double *mybuf)
buf_pos++;
}
//NP compute size
for(int ib = 0; ib < n_int_scalar_containers_; ib++)
{
//write owner data
fprintf(fp,"SCALARS %s int 1\nLOOKUP_TABLE default\n",int_scalar_container_names_[ib]);
m = buf_pos;
/*NL*/ //fprintf(screen,"size_one %d\n",size_one);
for (int i = 0; i < n; i++)
{
fprintf(fp,"%d\n",static_cast<int>(mybuf[m]));
m += size_one;
}
buf_pos++;
}
//NP compute size
for(int ib = 0; ib < n_scalar_containers_; ib++)
{

View File

@ -74,6 +74,9 @@ class DumpMeshVTK : public Dump {
class ScalarContainer<double> **min_active_edge_dist_;
// general implementation
class ScalarContainer<int> ***int_scalar_containers_;
char **int_scalar_container_names_;
int n_int_scalar_containers_;
class ScalarContainer<double> ***scalar_containers_;
char **scalar_container_names_;
int n_scalar_containers_;

13
src/error.cpp Normal file → Executable file
View File

@ -154,6 +154,19 @@ void Error::fix_error(const char *file, int line, Fix *fix,const char *str)
exit(1);
}
/* ----------------------------------------------------------------------
called by one proc in world
only write to screen if non-NULL on this proc since could be file
------------------------------------------------------------------------- */
void Error::fix_warning(const char *file, int line, Fix *fix, const char *str, int logflag)
{
if (screen) fprintf(screen,"WARNING: Fix %s (id %s): %s (%s:%d)\n",fix->style,fix->id,str,file,line);
if (logflag && logfile) fprintf(logfile,"WARNING: Fix %s (id %s): %s (%s:%d)\n",
fix->style,fix->id,str,file,line);
}
void Error::compute_error(const char *file, int line, Compute *compute,const char *str)
{
MPI_Barrier(world);

Some files were not shown because too many files have changed in this diff Show More