Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/test/python/test_bridging.py
bca23deb2955features/cadd

Last Change
Raw File

test_bridging.py
View Options

#!/usr/bin/python3
import pylibmultiscale as lm
import pytest
import numpy as np
################################################################
# geometrical properties
################################################################
r0 = 1.23157195294235
nb_atoms = 1000
LMD = nb_atoms*r0
el_size = 8.1*r0
nb_el = 300
LFE = el_size*nb_el
bridging_size = el_size*10
rcut = 2.2*r0
quality = 1e-3
epsilon = 1.6567944e-21/6.94769e-21 # kcal.mol^-1
sigma = 1.1 # angstrom
################################################################
lm.loadModules()
Dim = 3
lm.print_trace.fset(True)
comm = lm.Communicator.getCommunicator()
nb_proc = comm.getNumberFreeProcs()
comm.addGroup('md', nb_proc)
md_group = lm.Communicator.getGroup('md')
fe_group = md_group
md_geom = lm.Ball(Dim, 'md_geom')
md_geom.params.center = [0]
md_geom.params.radii = [0, LMD+bridging_size]
md_geom.init()
fe_geom = lm.Ball(Dim, 'fe_geom')
fe_geom.params.center = [0]
fe_geom.params.radii = [LMD, LMD+LFE]
fe_geom.init()
# coupling geom
bridge_geom = lm.Ball(Dim, 'bridge_geom')
bridge_geom.params.center = [0]
bridge_geom.params.radii = [LMD, LMD+bridging_size-el_size]
bridge_geom.init()
# coupling boundary geom
bound_geom = lm.Ball(Dim, 'bound_geom')
bound_geom.params.center = [0]
bound_geom.params.radii = [LMD+bridging_size-el_size, LMD+bridging_size]
bound_geom.init()
gManager = lm.GeometryManager.getManager()
gManager.addGeometry(md_geom)
gManager.addGeometry(fe_geom)
gManager.addGeometry(bridge_geom)
gManager.addGeometry(bound_geom)
lm.spatial_dimension.fset(1)
@pytest.fixture
def managers():
yield
lm.FilterManager.destroy()
lm.ActionManager.destroy()
lm.DomainMultiScale.destroy()
def test_create_domain(managers):
domMD = lm.DomainMD1D('md', md_group)
domMD.params.rcut = rcut
domMD.params.r0 = r0
domMD.params.mass = 39.95
domMD.params.timestep = 10
domMD.params.epsilon = 1.6567944e-21/6.94769e-21
domMD.params.sigma = 1.1
domMD.params.domain_geometry = "md_geom"
domMD.init()
lm.DomainMultiScale.getManager().addObject(domMD)
domFE = lm.DomainMeca1D('fe', fe_group)
domFE.params.elem_size = el_size
domFE.params.domain_geometry = "fe_geom"
domFE.params.mass = 39.95
domFE.params.r0 = r0
domFE.params.rcut = rcut
domFE.params.epsilon = 1.6567944e-21/6.94769e-21
domFE.params.sigma = 1.1
domFE.params.timestep = 10
domFE.init()
lm.DomainMultiScale.getManager().addObject(domFE)
return domMD, domFE
@pytest.fixture
def domains(managers):
return test_create_domain(managers)
def lj_p(r):
return 24*epsilon*sigma**6*(r**6 - 2*sigma**6)/r**13
def test_md_force(domains):
domMD, domFE = domains
atoms = domMD.getContainer()
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(atoms)
pos0 = extract.evalOutput().array().copy()
alpha = 1e-4
# alpha = 0
beta = 4.546
# beta = 0
linear_disp = pos0 * alpha + beta
set_field = lm.StimulationField('set')
set_field.params.field = lm.displacement
set_field.compute(dofs=atoms, field=linear_disp)
domMD.updateGradient()
extract.params.field = lm.force
extract.compute(atoms)
force = extract.evalOutput().array()
print(force[:20].flatten())
r = r0*(1. + alpha)
assert np.fabs(lj_p(r) + lj_p(2*r) - force[0]) < 1e-10
assert np.fabs(lj_p(r) + lj_p(2*r) + lj_p(-r) - force[1]) < 1e-10
assert (np.fabs(force[2:-2]) < 1e-10).any()
assert (np.fabs(force[1] + force[-2]) < 1e-10).any()
assert (np.fabs(force[0] + force[-1]) < 1e-10).any()
def test_md_force_pbc(domains):
domMD, domFE = domains
domMD.params.periodic = True
domMD.init()
atoms = domMD.getContainer()
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(atoms)
pos0 = extract.evalOutput().array().copy()
alpha = 1e-4
# alpha = 0
beta = 4.546
# beta = 0
linear_disp = pos0 * alpha + beta
set_field = lm.StimulationField('set')
set_field.params.field = lm.displacement
set_field.compute(dofs=atoms, field=linear_disp)
domMD.updateGradient()
extract.params.field = lm.force
extract.compute(atoms)
force = extract.evalOutput().array()
print(force[:20].flatten())
assert (np.fabs(force) < 1e-10).any()
def test_bridge_init(domains):
domMD, domFE = domains
assert domMD is not None
assert domFE is not None
bridge = lm.Bridging("test")
bridge.params.geometry = "bridge_geom"
bridge.init(domMD.getContainer(), domFE.getContainer())
return domMD, domFE, bridge
@pytest.fixture
def bridge_domain(domains):
return test_bridge_init(domains)
def test_bridge_constant_field_interpolation(bridge_domain):
domMD, domFE, bridge = bridge_domain
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(bridge.meshList)
pos0 = extract.evalOutput().array()
set_field = lm.StimulationField('set')
set_field.params.field = lm.position
set_field.compute(dofs=bridge.meshList, field=np.ones_like(pos0))
extract.params.field = lm.position
extract.compute(bridge.meshList)
assert (extract.evalOutput().array() == np.ones_like(pos0)).all()
out = lm.ContainerArrayReal(bridge.pointList.castContainer().size())
bridge.mesh2Point(lm.position, out)
assert (out.array() == np.ones_like(out.array())).all()
def test_bridge_linear_field_interpolation(bridge_domain):
domMD, domFE, bridge = bridge_domain
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(bridge.meshList)
pos0 = extract.evalOutput().array()
set_field = lm.StimulationField('set')
set_field.params.field = lm.position
alpha = 1.45
beta = 4.546
linear_field = pos0 * alpha + beta
set_field.compute(dofs=bridge.meshList, field=linear_field)
out = lm.ContainerArrayReal(bridge.pointList.castContainer().size())
bridge.mesh2Point(lm.position, out)
extract.params.field = lm.position0
extract.compute(bridge.pointList)
pos0 = extract.evalOutput().array()
assert out.array().shape == pos0.shape
assert (np.abs(out.array() - (pos0 * alpha + beta)) < 1e-12).all()
def test_xiao_init(domains):
domMD, domFE = domains
xiao = lm.Xiao("coupler")
xiao.params.geometry = "bridge_geom"
xiao.params.boundary = "bound_geom"
xiao.params.grid_division = [10]
xiao.params.quality = quality
xiao.init(domMD, domFE)
return xiao
@pytest.fixture
def xiao(domains):
return test_xiao_init(domains)
def test_xiao_checkweights(xiao):
xiao = xiao
bridge = xiao.bridging_zone
w_points = xiao.evalOutput('weight-point').evalOutput().array()
w_mesh = xiao.evalOutput('weight-mesh').evalOutput().array()
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(bridge.meshList)
pos0_mesh = extract.evalOutput().array().copy()
extract.compute(bridge.pointList)
pos0_points = extract.evalOutput().array().copy()
expected_w_mesh = np.zeros_like(pos0_mesh)
expected_w_points = np.zeros_like(pos0_points)
for i in range(0, pos0_mesh.shape[0]):
expected_w_mesh[i] = (np.abs(pos0_mesh[i, 0]) -
LMD)/(bridging_size - el_size)
for i in range(0, pos0_points.shape[0]):
expected_w_points[i] = 1. - (np.abs(pos0_points[i, 0]) -
LMD)/(bridging_size - el_size)
expected_w_points[expected_w_points == 0] = quality
expected_w_mesh[expected_w_mesh == 0] = quality
print("point weights")
print(w_points)
print("expected")
print(expected_w_points)
assert (expected_w_points == w_points).any()
print("mesh weights")
print(w_mesh)
print("expected")
print(expected_w_mesh)
assert (expected_w_mesh == w_mesh).any()
def test_xiao_rhs(xiao):
xiao = xiao
bridge = xiao.bridging_zone
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.position0
extract.compute(bridge.meshList)
pos0_mesh = extract.evalOutput().array().copy()
extract.compute(bridge.pointList)
pos0_points = extract.evalOutput().array().copy()
alpha1 = 1.45
beta1 = 4.546
vel_points = pos0_points * alpha1 + beta1
set_field = lm.StimulationField('set')
set_field.params.field = lm.velocity
set_field.compute(dofs=bridge.pointList, field=vel_points)
alpha2 = 2.54
beta2 = 5.098
vel_mesh = pos0_mesh * alpha2 + beta2
set_field.compute(dofs=bridge.meshList, field=vel_mesh)
xiao.buildRHS(lm.velocity)
rhs = xiao.evalOutput("rhs").array()
print(rhs)
expected_rhs = pos0_points*alpha2 + beta2 - vel_points
assert (rhs == expected_rhs).any()
def test_xiao_weighted_masses(xiao):
xiao = xiao
bridge = xiao.bridging_zone
xiao.buildConstraintMatrix()
A = xiao.evalOutput("A").array().flatten()
lambdas_points = xiao.evalOutput("lambdas-point").array().flatten()
lambdas_mesh = xiao.evalOutput("lambdas-mesh").array().flatten()
extract = lm.ComputeExtract('extractor')
extract.params.field = lm.mass
extract.compute(bridge.meshList)
mass_mesh = extract.evalOutput().array().copy()
extract.compute(bridge.pointList)
mass_points = extract.evalOutput().array().copy()
w_points = xiao.evalOutput('weight-point').evalOutput().array()
w_mesh = xiao.evalOutput('weight-mesh').evalOutput().array()
expected_lambdas_mesh = w_mesh * mass_mesh
expected_lambdas_points = w_points * mass_points
assert (expected_lambdas_mesh == lambdas_mesh).any()
assert (expected_lambdas_points == lambdas_points).any()
def test_xiao_constrained_matrix(xiao):
xiao = xiao
xiao.buildConstraintMatrix()
A = xiao.evalOutput("A").array().flatten()
lambdas_point = xiao.evalOutput("lambdas-point").array().flatten()
lambdas_mesh = xiao.evalOutput("lambdas-mesh").array().flatten()
smatrix = xiao.bridging_zone.smatrix.toarray()
print(smatrix.shape)
expected_A = np.einsum("Ii, I, I->i", smatrix,
smatrix.sum(axis=1), 1./lambdas_mesh)
print(expected_A.shape)
print(A.shape)
expected_A += 1./lambdas_point
print(A)
print(expected_A)
assert (expected_A == A).any()
c4science · Help