From e9e078501baa8d97a7376b550601cfd1ed9e8aaf Mon Sep 17 00:00:00 2001 From: Claudio Perez <50180406+claudioperez@users.noreply.github.com> Date: Sun, 22 Jun 2025 14:19:06 -0700 Subject: [PATCH 1/3] organize, use `getResidual` to get residual --- .gitignore | 2 + ddm.py | 148 +++++++++++++++++++++----------- main.py | 27 +++--- matmul | Bin 16704 -> 0 bytes matmul.c => matmul/matmul.c | 0 matmul_f.c => matmul/matmul_f.c | 0 model.py | 4 +- partitions.py | 14 +-- 8 files changed, 125 insertions(+), 70 deletions(-) create mode 100644 .gitignore delete mode 100755 matmul rename matmul.c => matmul/matmul.c (100%) rename matmul_f.c => matmul/matmul_f.c (100%) diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d646835 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.pyc +__pycache__/ diff --git a/ddm.py b/ddm.py index 5574e27..9768457 100644 --- a/ddm.py +++ b/ddm.py @@ -1,36 +1,35 @@ import numpy as np import pickle from multiprocessing import Process, Queue -import openseespy.opensees as ops -# ----------------------------------------------------------------------------- -# Domain – master process: spawns workers, assembles Schur system, shuts down -# ----------------------------------------------------------------------------- + class Domain: def __init__(self, partitions): """ + master process: spawns workers, assembles Schur system, shuts down + partitions: list of PartitionBuilder instances """ self.partitions = partitions self.workers = [] # list of Process - self.to_subdomain = [] # list of Queue (to worker) - self.from_subdomain = [] # list of Queue (from worker) + self._send = [] # list of Queue (to worker) + self._recv = [] # list of Queue (from worker) self.interface_guess = {} # {(node, dof): value} - self.interface_order = [] # [(node, dof), …] in global order + self.interface_order = [] # [(node, dof), ...] in global order # Spawn one Subdomain worker per partition for part in partitions: - to_q = Queue() - from_q = Queue() + q_send = Queue() + q_recv = Queue() p = Process( target=Domain._subdomain_worker, - args=(pickle.dumps(part), to_q, from_q) + args=(pickle.dumps(part), q_send, q_recv) ) p.start() self.workers.append(p) - self.to_subdomain.append(to_q) - self.from_subdomain.append(from_q) + self._send.append(q_send) + self._recv.append(q_recv) def step(self): """ @@ -43,20 +42,23 @@ def step(self): "dirichlet": self.interface_guess.copy(), "neumann": {} } - for q in self.to_subdomain: + for q in self._send: q.put({"cmd": "solve", "data": interface_data}) - return [q.get() for q in self.from_subdomain] + + return [q.get() for q in self._recv] + def schur_update(self): """ Assemble and solve the global Schur system. """ - for q in self.to_subdomain: + for q in self._send: q.put({"cmd": "schur"}) + S_total = None g_total = None - for q in self.from_subdomain: + for q in self._recv: S_i, g_i = q.get() if S_total is None: S_total = S_i.copy() @@ -76,80 +78,105 @@ def unpack_interface_vector(self, flat_u): } def shutdown(self): - for q in self.to_subdomain: + for q in self._send: q.put("TERMINATE") for p in self.workers: p.join() + @staticmethod - def _subdomain_worker(serialized_partition, to_q, from_q): - import pickle - from ddm import Subdomain - from partitions import LeftQuadPartition, RightQuadPartition + def _subdomain_worker(serialized_partition, q_send, q_recv): partition = pickle.loads(serialized_partition) sd = Subdomain(partition) while True: - msg = to_q.get() + msg = q_send.get() cmd = msg['cmd'] if isinstance(msg, dict) and 'cmd' in msg else msg if cmd == "schur": S_i, g_i = sd.get_schur_data() - from_q.put((S_i, g_i)) + q_recv.put((S_i, g_i)) + elif cmd in ("exit", "TERMINATE"): break + else: raise RuntimeError(f"Unknown message to subdomain worker: {msg}") -# ----------------------------------------------------------------------------- -# Subdomain – local model and Schur contribution -# ----------------------------------------------------------------------------- + class Subdomain: def __init__(self, partition): self.partition = partition + if not hasattr(self.partition, 'node_tags'): raise AttributeError("Partition object lacks node_tags") - ops.model("BasicBuilder", "-ndm", 2, "-ndf", 2) - partition.populate(ops) - ops.system("FullGeneral") - ops.numberer("Plain") - ops.constraints("Plain") - ops.integrator("LoadControl", 1.0) - ops.algorithm("Linear") - ops.analysis("Static") + import opensees.openseespy as _ops + model = _ops.Model("BasicBuilder", "-ndm", 2, "-ndf", 2) + self.model = model + + partition.populate(model) + + model.system("FullGeneral") + model.numberer("Plain") + model.constraints("Plain") + model.integrator("LoadControl", 1.0) + model.algorithm("Linear") + model.analysis("Static") dofs = partition.get_dof_partition() - self.interior_dofs = dofs["interior"] + self.interior_dofs = dofs["interior"] self.interface_dofs = dofs["interface"] + + def apply_interface_conditions(self, interface_data): for (node, dof), val in interface_data.get("dirichlet", {}).items(): - ops.sp(node, dof, val) + self.model.sp(node, dof, val) + for (node, dof), val in interface_data.get("neumann", {}).items(): - ops.load(node, dof, val) + self.model.load(node, dof, val) + def get_schur_data(self): print("Assembling tangent") - rc = ops.analyze(1) + self.model.system("FullGeneral") + self.model.numberer("Plain") + self.model.constraints("Plain") + self.model.test("FixedNumIter", 1) + self.model.integrator("LoadControl", 1.0) + self.model.algorithm("Linear") + self.model.analysis("Static") + + rc = self.model.analyze(1, operation="increment") if rc != 0: raise RuntimeError(f"OpenSees assemble failed (rc={rc})") - print("Assembly complete, fetching K") - N = ops.systemSize() - K = np.array(ops.printA("-ret"), dtype=float).reshape((N, N)) + +# N = self.model.systemSize() + K = self.model.getTangent() # Build global residual - R_global = np.zeros(N) - free_node_list = [n for n in self.partition.node_tags if n not in self.partition.fixed_nodes] - free_nodes = {n: i*2 for i, n in enumerate(free_node_list)} - for node, (Fx, Fy) in self.partition.get_nodal_loads().items(): - if node in free_nodes: - base = free_nodes[node] - R_global[base] = Fx - R_global[base+1] = Fy + R_global = self.model.getResidual() + +# free_node_list = [ +# n for n in self.partition.node_tags +# if n not in self.partition.fixed_nodes +# ] +# free_nodes = {n: i*2 for i, n in enumerate(free_node_list)} + +# for node, (Fx, Fy) in self.partition.get_nodal_loads().items(): +# if node in free_nodes: +# base = free_nodes[node] +# R_global[base] = Fx +# R_global[base+1] = Fy + + +# print(R_global) +# print(f"getResid = {self.model.getResidual()}") I, G = self.interior_dofs, self.interface_dofs + K_II = K[np.ix_(I, I)] K_IG = K[np.ix_(I, G)] K_GI = K[np.ix_(G, I)] @@ -158,7 +185,26 @@ def get_schur_data(self): R_G = R_global[G] KII_inv_KIG = np.linalg.solve(K_II, K_IG) - KII_inv_RI = np.linalg.solve(K_II, R_I) + KII_inv_RI = np.linalg.solve(K_II, R_I) S_local = K_GG - K_GI @ KII_inv_KIG - g_local = R_G - K_GI @ KII_inv_RI + g_local = R_G - K_GI @ KII_inv_RI return S_local, g_local + + +if __name__ == "__main__": + from partitions import LeftQuadPartition, RightQuadPartition + left = LeftQuadPartition() + right = RightQuadPartition() + + print("Starting domain setup") + domain = Domain([left, right]) + domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] + domain.interface_guess = {key: 0.0 for key in domain.interface_order} + + print("Calling schur_update") + u_gamma = domain.schur_update() + + print("schur_update completed, interface_guess:", domain.interface_guess) + domain.shutdown() + print("Domain shut down") + diff --git a/main.py b/main.py index 7cf829f..2b32f6a 100644 --- a/main.py +++ b/main.py @@ -1,15 +1,18 @@ from ddm import Domain from partitions import LeftQuadPartition, RightQuadPartition -left = LeftQuadPartition() -right = RightQuadPartition() -print("Starting domain setup") -domain = Domain([left, right]) -print("Domain initialized") -domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] -domain.interface_guess = {key: 0.0 for key in domain.interface_order} -print("Calling schur_update") -u_gamma = domain.schur_update() -print("schur_update completed, interface_guess:", domain.interface_guess) -domain.shutdown() -print("Domain shut down") +if __name__ == "__main__": + left = LeftQuadPartition() + right = RightQuadPartition() + + print("Starting domain setup") + domain = Domain([left, right]) + domain.interface_order = [(2,1), (2,2), (5,1), (5,2), (8,1), (8,2)] + domain.interface_guess = {key: 0.0 for key in domain.interface_order} + + print("Calling schur_update") + u_gamma = domain.schur_update() + + print("schur_update completed, interface_guess:", domain.interface_guess) + domain.shutdown() + print("Domain shut down") diff --git a/matmul b/matmul deleted file mode 100755 index ac29594c075f633ad637497b8f42f9de64091e79..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16704 zcmeHOeQ;dWb-ybiBZK9Yi50*k_~Bq;16oU#?6D1k^~W<$K{k#pkSXMsblk+#(J=CKHMk%gmnzhF z`#bmD^Y-mym8NO)SN3Z5oqK-gF8%m=DhS&_nyApl` zxMr)AU8(#^`c$V%cybPW=^S`s4*Va1YxsH7HUO2%rD>}YZp?xI9Jq#`H{A)KlAkQ_ zWnN|QQT!H(B|;0@LLP5oAZulWX=bcY-(zN9|6^t}6pItu*B?rYJrT=_B_gS4IJ8#` zSpDfpC~W?5DwcGZK0sDPgsfCdBtr3cst?}KS^LOk+A`#w(^jphNWHceJ+=ljRE_+jd@X)Hd8Iy0xjzNval1GcYoHnJs|tJh z&z?iqsuc$npIs80Pkz4#o)XWm9Op+>V-2eLlW|Sq^+9PS7-#TE`a60J1<^MY+IOM|V z-Z=fxg;N=v##}hZh-61y_%*ByI_|=+b>Syn_yQMx+J!H4;ni$40@VmqBT$V%H3HQL z{O%+0OaILu7{eFpjeOm|tP;Zb^0-w~JZ%iWR{yFjytrm3;NtRbT>VRfg7o`IGW}Mu zSd5NJo)#w4KP~aJaGIVd@w8x?eyPOMf@S)4(_bs`v_P5uQi-R9$#lBJ(}HBW zx5U%JYkH@{7nk3R+D83Lw+LC^;_`0=%K7I!{8v5v(;ogw5C6D_PkQ*Mhu`htA1vcr z|GhJJX19@h%NU-T+0m`Hj<>#HjC8z$uoauypn549@Gqs+F&b#3;~^l%NZn2%>u0Sd zG~LH3kEWYnob)dpBE44?h5kyZzh)Vcl~;`1jB)IpdyQkWwTADwapsD(7!J-i+YaiB zlXl$Lp4xY)<3or=_KqH7xMPB{M((_|z{qzT2UEyhDHaQ1ME`i*o8Wzqz}6|B?5CfA zkTN~6LuYHr?LU{-&-ssx8wZVJleI?fMB!yPF>?A0@DliK8=O6%U!Z$>yDhI@@E?&H zBMV+92msVa7E8v+>$8UMw2|+bF?>HVa*f0K1s^o-($DRCorG_^>eM2)zmUsLozQ2c z)QvL5Pd~MPgXLj3O@IrBCg9ZCIa0 zewOk_Mn@kh6*hc9Gs+5_{ys|MGPnFDWg2-AHT=nOH2i0=jz-9i(x@wReON5ob(k71 z_}1%}_r3VPbXV_o_u=JY@fBLv7hXlKbVL3F+tG%jqh%31iI|?ui=1A_*G=7l#%Op8 z*KhpiuX6f@yk0;{U63tW7}-J1;6E~v(=X-qnVf!Z`e!l#lemiMmvM>u-yKIhb+q*q zfp16AE1xH)J?Bc1M-=isQzafPK{MAAQqP4;y~#2?wAMq%pucVQqW=$&K~cl{lrLX* zJB-(($)1IywG|62ccjAswfiJ&-?^P^C;bOH5S?ND9B%X6UHT>efm+1mgSNrv9D~nOZAR*HO6&b!kyeYXwcgRXUux;6|A9hH#b1Hm3H>m! zZ{Pn4{o5*BYu~8d_j?ZK8uR+$vKiy3F>-ee;(yY}9Yc|)$Y*Xsx-9%F{0-}e(Q}TX z3#``<@B1ruzv82QXetvn?X0Hw3G(l{wiK&j!7RN5^^ z6$Y<8s*q{a>C*6&5@X(aasqwe(-`HL!X|RZQMbhqQrD7`*jHxJsVT@PX*O)!u~-fe ztFEvc0(;7lYZPZyLb}ydjX*U5)d*B0P>nz}0>5np=uHfz{8Qr7Zoy;bBBWW!a)$G?&nR~E06eQG^hqXjDV`;q+EreX$iVeggqzN0B0rajB z+ia(;8#uSx$+c||V!KnywjG^XI+T1Yvq5W_UxHQP8{lldJ{alCTCr48qnBnImIvD6 zQ7secql(F--ogZ*4f<}Xy^G#)VgV@Dy<03k16lwb1Dypu33_m*SR4R->b+ud67=ee z#Uj0Qu|PG@&R-UbLD0u96^oD3!Uc)m=~D82VW;p7Hu`Q_aCQACWQ5a8SPedJBZVzj z%hx)$BHNDN58tuN77H7>tLj1$A@34n#~RKYWAZ>c3o&4TZBAOCJ#IEQ^31$?{AjLPdV}nzt@OOdk$2Qe2Y{-n8i8sAsu8G0pc;W{1ga6JMxYvj-?s?x zcWL|`n%YLlZ{H{u3pgNZ(T~?mGX41tHoZ~t{Jq;vipLUPrX`BqNhKxO9HusEyayp_77@M{%S)e%jKXQIr+fcn16-d*B*Wu-v-7fSrC8GVmR z>5u{#KdCCh^*(E>$^9E;x%`_I&U$#JE$#M|jMXwg^zkQG#PyyDI!pA@>=~MgG zh(lMZ^Z3N|!rOmNNRy8_-_f%t`i|L$QNrVp z@fl<)eaCzi0#XUT1-Qn>F8Z9rF&hP4AhI36>AU9^h0{ATN}GUJ^1lmsqqtEt2Bn&F z`9spLStNMgW5J_y_<3;-{7>e9klg~gR%;-IujWI&G1Cg|5dbn)HX04|iP9JLrj;=J z=o@=FBoI!Sd*Z3yP}~e#sdUB+Wd}uHDuIWqRwNu~sjNZY`NzyqIvv_;Mv_)~uZX5Y ziHI4_CK7vL;^ItbT4kzuSb%gn>Lv;610{NO9`m*PS$=hs^TegkvDx zMZzH~B;@^~V&x|)9L@H!GZhYO*M}U~2N-acf*y9S7Eh$ok42oL6CVBY0}p06lQR25 z$uJ#;=zIW=Q=k?x|Dc+wK*ZyBICCAo*$(59qz>QBxAHL zbKCQ}kg3LsihAuo2mF1k*}49_uVK1E*{S9KN9rGfoc1(q&-)&x2PwEnR7cA4{qz1K z28@=nY|r~Bro3;0icHEKv*Po}(Vm3uc^|@bnbPF?vmVnYA*VeK^SsYt8dUb)`YTwB zAVV=_d){|3<$V{I@4f%$mHj5=@DLS*#FW;Xl*pc9=C%J0Fp3HLcmB_+|F^uqg^Enx z_J0AW+g?+4OvjWR6D4+v`GLp&v{GQo`g7Tne|Pzl%7E!1HdK`TG5uSQJ+IH0-m2os zZrG0bS3LH-US`VA4_MzDKl*T%>d);jG^MF%L9tA@{h0qTGH&~VW=oiMv!bHj@_zx` zZ6CbFR%Ci730=us{&~n2;>Y`SelEi2eZBh3Gkp)XewRJ(0|%P{BC(#$m3YQKgn;@k zm(Tk%KEKW9w7D7Cj@$9CVN18j_Pk#kSOpMi)F#fA?U~MlvD;n*m1>Jkxg2)DcFfZ| zTf(?}Eoe&`mAxXYE;J`+XBMgQy?Shq4{FO2)s>Q)I(^R;$!kzKT%n#zeG({bQMmp* l|IP=-b(e Date: Sun, 22 Jun 2025 14:20:08 -0700 Subject: [PATCH 2/3] update readme --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8885b8b..863753b 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ An example mesh showing quadrilateral elements and node numbering: - Python - Parallel Domain Decomposition -- Finite Element concepts +- Finite Element concepts with [xara](https://xara.so) ## File Structure From 96f347e708f82b5e733f42e90bde7471211578b1 Mon Sep 17 00:00:00 2001 From: Claudio Perez <50180406+claudioperez@users.noreply.github.com> Date: Sun, 22 Jun 2025 14:23:31 -0700 Subject: [PATCH 3/3] revert organizing --- matmul/matmul.c => matmul.c | 0 matmul/matmul_f.c => matmul_f.c | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename matmul/matmul.c => matmul.c (100%) rename matmul/matmul_f.c => matmul_f.c (100%) diff --git a/matmul/matmul.c b/matmul.c similarity index 100% rename from matmul/matmul.c rename to matmul.c diff --git a/matmul/matmul_f.c b/matmul_f.c similarity index 100% rename from matmul/matmul_f.c rename to matmul_f.c