Skip to content

Commit

Permalink
Merge pull request #7 from Argonne-National-Laboratory/network_fix
Browse files Browse the repository at this point in the history
Fix issues with the spring based "system" model
  • Loading branch information
reverendbedford authored Mar 18, 2021
2 parents 9fa4929 + 3d55992 commit ac31291
Show file tree
Hide file tree
Showing 8 changed files with 595 additions and 21 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -9,5 +9,5 @@ jobs:
- run: pip3 install --user wheel
- run: pip3 install --user -r requirements.txt
- run: nosetests3
- run: pylint3 srlife
- run: pylint srlife

45 changes: 29 additions & 16 deletions srlife/spring.py
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,8 @@ def __init__(self, k):
k spring stiffness
"""
self.k = k
self.state_np1 = None
self.state_n = None

def force_and_stiffness(self, i, d):
"""
Expand Down Expand Up @@ -238,7 +240,7 @@ def remove_rigid(self):
Remove rigid links by merging the nodes in place
"""
while True:
for i,j, edge in self.edges(data=True):
for i,j,key, edge in self.edges(data=True,keys=True):
# Can't handle duplicate springs
if edge['object'] == "rigid":
if self.number_of_edges(i,j) != 1:
Expand All @@ -247,10 +249,15 @@ def remove_rigid(self):
bcj = 'bc' in self.nodes[j]
if bci and bcj:
raise RuntimeError("Cannot merge two nodes with BCs!")
self.remove_edge(i,j)
if bci or ((not bci) and (not bcj)):
self.remove_edge(i,j,key=key)

if i < j:
if bcj:
raise RuntimeError("Internal error: deleting BC")
self.replace(nx.contracted_nodes(self, i, j))
else:
if bci:
raise RuntimeError("Internal error: deleting BC")
self.replace(nx.contracted_nodes(self, j, i))
break
else:
Expand Down Expand Up @@ -295,6 +302,7 @@ def validate_solve(self):
Conditions:
1) All dummy edges removed
2) At least one fixed condition
3) All nodes sequentially numbers from zero
"""
# Check that we have times
if self.times is None:
Expand Down Expand Up @@ -347,7 +355,7 @@ def solve(self, i, nthreads = 1):
d = newton(lambda x: self.RJ(x, nthreads = nthreads),
self.displacements[self.dmap[self.free]], verbose = self.verbose,
rel_tol = self.rtol, abs_tol = self.atol, miters = self.miter)

# Store the displacements, for fun
self.displacements[self.dmap[self.free]] = d
self.displacements[self.dmap[self.fixed]] = self.fixed_displacements
Expand All @@ -372,6 +380,7 @@ def RJ(self, d, nthreads = 1):
res = list(p.map(lambda e: self.fj(dall, *e), self.edges(data=True)))
else:
res = list(map(lambda e: self.fj(dall, *e), self.edges(data=True)))

Fint = sum(r[0] for r in res)
J = sum(r[1] for r in res)

Expand All @@ -385,19 +394,21 @@ def fj(self,dall,i,j,edge):
"""
Calculate the force and jacobian contributions for a single edge
"""
ii = self.dmap[i]
jj = self.dmap[j]
ss = np.sign(ii-jj)
Fint = np.zeros((len(self.nodes),))
J = np.zeros((len(self.nodes),len(self.nodes)))
f, k = edge['object'].force_and_stiffness(self.i, dall[self.dmap[j]]
- dall[self.dmap[i]])
Fint[self.dmap[i]] += -f
Fint[self.dmap[j]] += f

J[self.dmap[i],self.dmap[i]] += k
J[self.dmap[i],self.dmap[j]] += -k
J[self.dmap[j],self.dmap[i]] += -k
J[self.dmap[j],self.dmap[j]] += k
f, k = edge['object'].force_and_stiffness(self.i, (dall[jj] - dall[ii])*ss)
Fint[ii] += -f
Fint[jj] += f

J[ii,ii] += k
J[ii,jj] += -k
J[jj,ii] += -k
J[jj,jj] += k

return Fint, J, edge['object'].state_np1
return Fint*ss, J, edge['object'].state_np1

def dof_maps(self, i):
"""
Expand All @@ -411,7 +422,9 @@ def dof_maps(self, i):
forces = []
displacements = []
time = self.times[i]
for node in self.nodes:
nodes = []
for node in sorted(self.nodes):
nodes.append(node)
if 'bc' in self.nodes[node]:
if self.nodes[node]['bc'][0] == 'displacement':
fixed.append(node)
Expand All @@ -425,7 +438,7 @@ def dof_maps(self, i):
free.append(node)
forces.append(0)

nodes = np.array(self.nodes, dtype=int)
nodes = np.array(nodes, dtype=int)
total = (np.array(range(0, max(nodes)+1), dtype = int)*0)-1
rnums = np.array(range(len(nodes)), dtype = int)
total[nodes] = rnums
Expand Down
6 changes: 3 additions & 3 deletions srlife/structural.py
Original file line number Diff line number Diff line change
Expand Up @@ -266,7 +266,7 @@ def solve(self, tube, i, state_n, dtop):
dtop, self.solver_options)
else:
raise ValueError("Unknown dimension %i" % tube.ndim)

return state_np1

def init_state(self, tube, mat, i = None):
Expand Down Expand Up @@ -694,7 +694,7 @@ def solve(self, t_n, t_np1, p):
# Calculate the initial residual and jacobian
R = self.residual(p)
nR0 = la.norm(R[self.kdofs])

# Printing, if you want
if self.options['verbose']:
print("Iter\tnR\t\tnR/nR0")
Expand Down Expand Up @@ -990,7 +990,7 @@ def calculate_stress_update(self):
).transpose(2,0,1)
self.state_np1.tangent = A_np1.reshape((self.state_np1.ne,
self.state_np1.nqi, 3,3,3,3)).transpose(2,3,4,5,0,1)

def get_mat(thing):
"""
Small helper to wrap NEML for pickling issues
Expand Down
3 changes: 3 additions & 0 deletions srlife/system.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,9 @@ def make_network(self, model, smat, ssolver):
Setup the complete spring network, given the receiver and
material objects
This function must produce a network where the nodes are numbered
topologically where bottom_tube < top_tube < panel < start
Parameters:
model fully-defined receiver object
smat structural material model
Expand Down
Empty file added test/solvers/spring/__init__.py
Empty file.
174 changes: 174 additions & 0 deletions test/solvers/spring/altsolver.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
import numpy as np
import numpy.linalg as la
import networkx as nx

class Network(nx.MultiGraph):
"""
Simplified version of panel network where everything is an elastic spring
Boundary conditions are always fixed on the bottom and free on the top
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
if len(args) == 0:
self.add_node("top", distance = 0)
self.cpanel = 0

def add_panel(self, k):
"""
Add a panel to the network with stiffness k
Parameters:
k: either a float stiffness or "rigid" or "disconnect"
Returns:
id of the new panel
"""
name = self._panel_name(self.cpanel)
self.cpanel += 1
self.add_node(name, distance = 1)
self.add_edge("top", name, stiffness = k)

return self.cpanel - 1

def add_tube(self, panel, kconnect, k):
"""
Add a tube to a panel
Parameters:
panel: panel id
kconnect: either a float stiffness or "rigid" or "disconnect"
k: float stiffness
"""
if panel >= self.cpanel:
raise ValueError("Panel %i does not exist!" % panel)

name = self._tube_name(panel, self._next_tube(panel))
self.add_node(name+"_top", distance = 2)
self.add_node(name+"_bottom", distance = 3)

self.add_edge(self._panel_name(panel), name+"_top", stiffness = kconnect)
self.add_edge(name+"_top",name+"_bottom", stiffness = k)

def simplify_network(self):
"""
Simplify the tube network and return a list of subproblems
"""
# Remove all rigid nodes
self._remove_rigid()

# Split at all disconnects
return self._split_disconnect()

def solve_network(self, tforce):
"""
Impose the standard BCs and solve for the displacements at each node
Parameters:
tforce: thermal force on the bars
"""
n = len(self)
K = np.zeros((n,n))
d = np.zeros((n,))
f = np.zeros((n,))

# Make the dof maps
self.backmap = {}
for i,name in enumerate(self.nodes()):
self.backmap[name] = i

# Fix all the ones with "_bottom" in the name to zero
fixed = []
for i,name in enumerate(self.nodes()):
if name[:4] == "tube" and name.split('_')[1] == "bottom":
fixed.append(i)
free = sorted(list(set(list(range(n))) - set(fixed)))

# Impose the thermal force and the stiffness matrix
for a,b,data in self.edges(data=True):
dofs = [self.backmap[a],self.backmap[b]]
s = np.sign(self.nodes[b]['distance']-self.nodes[a]['distance'])
K[np.ix_(dofs,dofs)] += (
data['stiffness']*np.array([[1,-1],[-1,1]]))

if a[:4] == 'tube' and b[:4] == 'tube':
f[dofs] += tforce * np.array([1,-1]) * s

# Solve
fprime = f[free] - np.dot(K[np.ix_(free,fixed)], d[fixed])
Kprime = K[np.ix_(free,free)]
dprime = la.solve(Kprime,fprime)
d[free] = dprime

# Backsub the displacements
for i, (name,data) in enumerate(self.nodes(data=True)):
data['displacement'] = d[i]

def _remove_rigid(self):
"""
Remove all rigid connections by deleting the edge and joining the
nodes
"""
while True:
for a,b,key,data in self.edges(data=True,keys=True):
if data['stiffness'] == "rigid":
self.remove_edge(a,b,key=key)
# Need to keep the one marked "tube"
if a[:4] == "tube":
ng = nx.contracted_nodes(self, a, b)
else:
ng = nx.contracted_nodes(self, b, a)
self.clear()
self.update(ng)
break
else:
break

def _split_disconnect(self):
"""
Split the graph at all disconnects
"""
# First just split
while True:
for i,j,key,data in self.edges(keys = True, data = True):
if data['stiffness'] == 'disconnect':
self.remove_edge(i,j,key = key)
break
else:
break

# Setup the new structures
base = [Network(self.subgraph(c)) for c in nx.connected_components(self) if len(c) > 1]
nstructs = []
for b in base:
ns = [n[:4] for n in b.nodes()]
if "tube" in ns:
nstructs.append(b)

return nstructs

def _panel_name(self, k):
"""
Turn an integer into a panel id
Parameters:
k: panel integer id
"""
return "panel-%i" % k

def _tube_name(self, panel, tube):
"""
Turn integer panel and tube ids into a base tube string id
"""
return "tube-%i-%i" % (panel, tube)

def _next_tube(self, panel):
"""
Figure out what the next tube id is for a panel
"""
start = 0
for n in self.neighbors(self._panel_name(panel)):
if n != "top":
start += 1

return start
Loading

0 comments on commit ac31291

Please sign in to comment.