sph_benchmark_dropletΒΆ

../_images/sph_benchmark_droplet.png
#!/usr/bin/env python-mpacts
from mpacts.commands.continuum.sph_eos import EOS_Monaghan
from mpacts.commands.monitors.progress import ProgressIndicator
from mpacts.commands.onarrays.setvalue import SetValueCommand
from mpacts.commands.time_evolution.integration import LeapFrog_Translation
from mpacts.commands.time_evolution.integration import ForwardEuler_Generic
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.models.continuum.sph import SPH_2D
from mpacts.core.arrays import ArrayManager
from mpacts.core.arrays import create_array
from mpacts.core.geometrytags import GeometryTag_SPH
from mpacts.core.simulation import Simulation
from mpacts.geometrygenerators.boolean3Dshapegenerator import Boolean3DShapeGenerator
from mpacts.io.vtk_writer import VTKSerialWriter
import mpacts.core.command as cmd
import random
import numpy as np

### Evolution of an elliptical drop ###
# Results agree with "Simulating Free Surface Flows with SPH" (Monaghan 1992):
# b (Semi-Major Axis) of the elliptical drop (Table I):
# Time [s]      b in Theory       b in this simulation
# 0.0008        1.083             1.0809
# 0.0038        1.44              1.4302
# 0.0076        1.95              1.9376


def create_sph_droplet(name):
    mgr = ArrayManager(name)
    x = create_array("Point", "x", mgr)
    v = create_array("Vector", "v", mgr)
    h = create_array("Scalar", "h", mgr)
    m = create_array("Scalar", "m", mgr)
    rho = create_array("Scalar", "rho", mgr)
    drho = create_array("Scalar", "drho", mgr)
    rhoNorm = create_array("Scalar", "rhoNorm", mgr)
    rhoSum = create_array("Scalar", "rhoSum", mgr)
    p = create_array("Scalar", "p", mgr)
    F = create_array("Vector", "F", mgr)
    mgr.add_child(GeometryTag_SPH())
    vxt = 100.0
    vyt = 100.0

    positions = Boolean3DShapeGenerator()
    positions.addVerticesCirc( center=( 0,0,0)
                             , r=1.00
                             , plane='xy'
                             , dist=0.04
                             , boolean='None' )

    xp = positions.getVertices()
    x.set_array(xp)
    vp = np.array(positions.getVertices())
    vp[:, 0] *= vxt * (-1)
    vp[:, 1] *= vyt
    v.set_array( list( map(tuple, vp) ) )
    mgr["F"].set_array((0, 0, 0))
    mgr["rho"].set_array(1000.0)
    mgr["drho"].set_array(0)
    mgr["rhoNorm"].set_array(0)
    mgr["rhoSum"].set_array(0)
    m.set_array(1.60)
    h.set_array(0.060)
    return mgr

mysim = Simulation("simulation", timestep=1e-6)
# print mysim.print_tree()
r = 0.1
B = 2.8e8
pc1 = create_sph_droplet("drop")
mysim.add_child(pc1)
SetValueCommand("zero_forces_pc1", mysim, value=(0, 0, 0), array=pc1["F"])
SetValueCommand("zero_rho_pc1", mysim, value=0, array=pc1["drho"])
EOS_Monaghan("eos_pc1", mysim, B=B, rho0=1000.0, pc=pc1)
LeapFrog_Translation("TranslationTigr", mysim, pc=pc1)
ForwardEuler_Generic("rho_pc1", mysim, pc=pc1, x=pc1["rho"], dx=pc1["drho"])

cmodel = SPH_2D("sph1", pc1=pc1, pc2=pc1, viscosity=1e5)
noc = create_array("unsigned", "number_of_contacts", pc1)
noec = create_array("unsigned", "number_of_effective_contacts", pc1)
SetValueCommand("zero_noc",mysim,value=0,array=pc1["number_of_contacts"])
SetValueCommand("zero_enoc",mysim,value=0,array=pc1["number_of_effective_contacts"])

cd = MultiGridContactDetector( "sph2", mysim, cmodel=cmodel, update_every=40, keep_distance = 0.04 )

wrt = VTKSerialWriter( "test",mysim("loop_cmds/output_cmds"), data=mysim, filename="sph"
                     , gate=cmd.ExecuteEveryNTimes(N=250)
                     , select_all = True )

ProgressIndicator("PrintProgress", mysim, print_interval=1)

mysim.run_until(0.01)

#-----------------------------------------------------------------------------------------------------
#Verification that the output is identical to a reference value
from mpacts.tools.arrayhash import check_and_document_hash

check_and_document_hash( pc1['x'], '97e0e2ae3658b6a12f8d63a5c75c5b77', tolerance=1e-7 )
#-----------------------------------------------------------------------------------------------------