Skip to content

Commit c51feb5

Browse files
matthiasneunermatthiasneuner
committed
Merge branch 'matthias/formatting' into matthias/refactoring_rkpm_weak_dirichlet_locations
2 parents a03068d + a72c992 commit c51feb5

File tree

3 files changed

+705
-33
lines changed

3 files changed

+705
-33
lines changed

‎edelweissmpm/solvers/nqs.py‎

Lines changed: 195 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,10 @@
2424
# The full text of the license can be found in the file LICENSE.md at
2525
# the top level directory of EdelweissMPM.
2626
# ---------------------------------------------------------------------
27-
27+
from collections.abc import Callable
2828

2929
import edelweissfe.utils.performancetiming as performancetiming
30+
import numpy as np
3031
from edelweissfe.constraints.base.constraintbase import ConstraintBase
3132
from edelweissfe.journal.journal import Journal
3233
from edelweissfe.numerics.dofmanager import DofManager, DofVector
@@ -88,6 +89,10 @@ class NonlinearQuasistaticSolver(NonlinearImplicitSolverBase):
8889
"spec. absolute field correction tolerances": dict(),
8990
"failed increment cutback factor": 0.25,
9091
"fall back to quasi Newton after allowed residual growths": False,
92+
"line search": False,
93+
"line search after n iterations": 5,
94+
"line search every n iterations": 2,
95+
"line search alphas": [0.8, 1.0, 1.2],
9196
}
9297

9398
def __init__(self, journal: Journal):
@@ -490,41 +495,30 @@ def _newtonSolve(
490495
quasi_Newton = False
491496

492497
while True:
493-
PInt[:] = K_VIJ[:] = F[:] = PExt[:] = 0.0
494-
495-
self._prepareMaterialPoints(materialPoints, timeStep.totalTime, timeStep.timeIncrement)
496-
self._interpolateFieldsToMaterialPoints(activeCells, dU)
497-
self._interpolateFieldsToMaterialPoints(cellElements, dU)
498-
self._computeMaterialPoints(materialPoints, timeStep.totalTime, timeStep.timeIncrement)
499-
500-
self._computeCells(
501-
activeCells, dU, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
502-
)
503-
504-
self._computeElements(
505-
elements, dU, Un, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
506-
)
507-
508-
self._computeCellElements(
509-
cellElements, dU, Un, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
510-
)
511498

512-
self._computeParticles(
513-
particles, dU, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
499+
Rhs, K_VIJ, F, PInt, PExt = self._computeSystem(
500+
dirichlets,
501+
bodyLoads,
502+
distributedLoads,
503+
particleDistributedLoads,
504+
reducedNodeSets,
505+
elements,
506+
Un,
507+
activeCells,
508+
cellElements,
509+
materialPoints,
510+
particles,
511+
constraints,
512+
timeStep,
513+
theDofManager,
514+
dU,
515+
Rhs,
516+
K_VIJ,
517+
PInt,
518+
PExt,
519+
F,
514520
)
515521

516-
self._computeConstraints(constraints, dU, PInt, K_VIJ, timeStep)
517-
518-
PExt, K = self._computeBodyLoads(bodyLoads, PExt, K_VIJ, timeStep, theDofManager, activeCells)
519-
PExt, K = self._computeCellDistributedLoads(distributedLoads, PExt, K_VIJ, timeStep, theDofManager)
520-
521-
PExt, K = self._computeParticleDistributedLoads(
522-
particleDistributedLoads, PExt, K_VIJ, timeStep, theDofManager
523-
)
524-
525-
Rhs[:] = -PInt
526-
Rhs -= PExt
527-
528522
if iterationCounter == 0 and dirichlets:
529523
Rhs = self._applyDirichlet(timeStep, Rhs, dirichlets, reducedNodeSets, theDofManager)
530524
else:
@@ -587,9 +581,177 @@ def _newtonSolve(
587581
K_CSR = self._applyDirichletKCsr(K_CSR, dirichlets, theDofManager, reducedNodeSets)
588582

589583
ddU = self._linearSolve(K_CSR, Rhs, linearSolver)
584+
585+
if (
586+
iterationOptions["line search"]
587+
and iterationCounter > iterationOptions["line search after n iterations"]
588+
and iterationCounter % iterationOptions["line search every n iterations"] == 0
589+
):
590+
591+
alphas = iterationOptions["line search alphas"]
592+
593+
def _computeResidualCallback(ddU_linesearch: DofVector) -> DofVector:
594+
595+
dU_linesearch = dU.copy()
596+
dU_linesearch += ddU_linesearch
597+
598+
Rhs_, *rest = self._computeSystem(
599+
dirichlets,
600+
bodyLoads,
601+
distributedLoads,
602+
particleDistributedLoads,
603+
reducedNodeSets,
604+
elements,
605+
Un,
606+
activeCells,
607+
cellElements,
608+
materialPoints,
609+
particles,
610+
constraints,
611+
timeStep,
612+
theDofManager,
613+
dU_linesearch,
614+
Rhs,
615+
K_VIJ,
616+
PInt,
617+
PExt,
618+
F,
619+
)
620+
for dirichlet in dirichlets:
621+
Rhs_[self._findDirichletIndices(theDofManager, dirichlet, reducedNodeSets[dirichlet.nSet])] = (
622+
0.0
623+
)
624+
625+
return Rhs_
626+
627+
ddU = self._quadraticLineSearch(_computeResidualCallback, ddU, alphas)
628+
629+
else:
630+
if initialGuess is not None:
631+
# We only do one iteration with the initial guess, then we switch to standard Newton
632+
quasi_Newton = False
633+
590634
dU += ddU
591635
iterationCounter += 1
592636

593637
iterationHistory = {"iterations": iterationCounter, "incrementResidualHistory": incrementResidualHistory}
594638

595639
return dU, PInt, iterationHistory, newtonCache
640+
641+
@performancetiming.timeit("compute system")
642+
def _computeSystem(
643+
self,
644+
dirichlets: list[DirichletBase],
645+
bodyLoads: list,
646+
distributedLoads: list,
647+
particleDistributedLoads: list,
648+
reducedNodeSets: list,
649+
elements: list,
650+
Un: DofVector,
651+
activeCells: list,
652+
cellElements: list,
653+
materialPoints: list,
654+
particles: list,
655+
constraints: list,
656+
timeStep: TimeStep,
657+
theDofManager: DofManager,
658+
dU,
659+
Rhs,
660+
K_VIJ,
661+
PInt,
662+
PExt,
663+
F,
664+
):
665+
"""Compute the global residual vector and the global stiffness matrix."""
666+
667+
PInt[:] = K_VIJ[:] = F[:] = PExt[:] = 0.0
668+
669+
self._prepareMaterialPoints(materialPoints, timeStep.totalTime, timeStep.timeIncrement)
670+
self._interpolateFieldsToMaterialPoints(activeCells, dU)
671+
self._interpolateFieldsToMaterialPoints(cellElements, dU)
672+
self._computeMaterialPoints(materialPoints, timeStep.totalTime, timeStep.timeIncrement)
673+
674+
self._computeCells(activeCells, dU, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager)
675+
676+
self._computeElements(
677+
elements, dU, Un, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
678+
)
679+
680+
self._computeCellElements(
681+
cellElements, dU, Un, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager
682+
)
683+
684+
self._computeParticles(particles, dU, PInt, F, K_VIJ, timeStep.totalTime, timeStep.timeIncrement, theDofManager)
685+
686+
self._computeConstraints(constraints, dU, PInt, K_VIJ, timeStep)
687+
688+
PExt, K = self._computeBodyLoads(bodyLoads, PExt, K_VIJ, timeStep, theDofManager, activeCells)
689+
PExt, K = self._computeCellDistributedLoads(distributedLoads, PExt, K_VIJ, timeStep, theDofManager)
690+
691+
PExt, K = self._computeParticleDistributedLoads(particleDistributedLoads, PExt, K_VIJ, timeStep, theDofManager)
692+
693+
Rhs[:] = -PInt
694+
Rhs -= PExt
695+
696+
return Rhs, K, F, PInt, PExt
697+
698+
@performancetiming.timeit("line search")
699+
def _quadraticLineSearch(
700+
self,
701+
computeResidual: Callable[[DofVector], DofVector],
702+
ddU: DofVector,
703+
alphas: list,
704+
):
705+
"""Perform a quadratic line search to determine an optimal step length.
706+
707+
Parameters
708+
----------
709+
computeSystem
710+
A callback function to compute the system residual for a given solution increment.
711+
ddU
712+
The current solution increment vector.
713+
alphas
714+
The list of alpha values to be tried.
715+
716+
Returns
717+
-------
718+
DofVector
719+
The updated solution increment vector after line search.
720+
"""
721+
722+
R_trial_values = []
723+
ddU_linesearch = ddU.copy()
724+
725+
for alpha in alphas:
726+
ddU_linesearch[:] = ddU * alpha
727+
728+
Rhs_trial = computeResidual(ddU_linesearch)
729+
730+
R_trial = np.linalg.norm(Rhs_trial, np.inf)
731+
R_trial_values.append(R_trial)
732+
self.journal.message(
733+
" line search try with alpha {:7.4f}, ||R||∞ {:11.9e}".format(alpha, R_trial),
734+
self.identification,
735+
level=2,
736+
)
737+
738+
coeffs = np.polyfit(alphas, R_trial_values, 2)
739+
alpha_opt = -coeffs[1] / (2 * coeffs[0])
740+
741+
alpha = max(0.1, min(1.5, alpha_opt))
742+
self.journal.message(
743+
" line search selected alpha {:7.4f} from quadratic fit".format(alpha),
744+
self.identification,
745+
level=2,
746+
)
747+
748+
ddU_final = ddU * alpha
749+
Rhs_final = computeResidual(ddU_final)
750+
751+
self.journal.message(
752+
" line search final ||R||∞ {:11.9e}".format(np.linalg.norm(Rhs_final, np.inf)),
753+
self.identification,
754+
level=2,
755+
)
756+
757+
return ddU_final

‎examples/123_quadratic_line_search_test/gold.csv‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
-6.853366687659331880e-04
2+
-1.020726553815671345e+00
3+
0.000000000000000000e+00
4+
2.476748344325084149e-04
5+
-1.021164846196251030e+00
6+
0.000000000000000000e+00
7+
-1.042038827845916072e-04
8+
-1.021477845836610943e+00
9+
0.000000000000000000e+00
10+
-1.078132833985714347e-04
11+
-1.021633493290713846e+00
12+
0.000000000000000000e+00
13+
-1.898047682925268077e-04
14+
-1.021832338479531366e+00
15+
0.000000000000000000e+00
16+
-1.605551375300816124e-04
17+
-1.022031086970464742e+00
18+
0.000000000000000000e+00
19+
-9.463992530927549005e-05
20+
-1.022234850132560657e+00
21+
0.000000000000000000e+00
22+
-3.847418423268506797e-05
23+
-1.022309239675273318e+00
24+
0.000000000000000000e+00
25+
-2.531837157376833460e-04
26+
-1.021471584907200691e+00
27+
0.000000000000000000e+00
28+
1.824048464224312734e-04
29+
-1.021302149434530904e+00
30+
0.000000000000000000e+00
31+
-3.181616296050254961e-05
32+
-1.021455640587263902e+00
33+
0.000000000000000000e+00
34+
-3.551388313610596254e-05
35+
-1.021633100597821642e+00
36+
0.000000000000000000e+00
37+
-7.657605133675757854e-05
38+
-1.021837513099376649e+00
39+
0.000000000000000000e+00
40+
-5.741699162678655740e-05
41+
-1.022083073545886345e+00
42+
0.000000000000000000e+00
43+
-3.051019645323016624e-05
44+
-1.022284054407304854e+00
45+
0.000000000000000000e+00
46+
-1.160920731194348766e-05
47+
-1.022354263954230236e+00
48+
0.000000000000000000e+00
49+
2.531848389800057421e-04
50+
-1.021471585047997399e+00
51+
0.000000000000000000e+00
52+
-1.824049725590419242e-04
53+
-1.021302149338856324e+00
54+
0.000000000000000000e+00
55+
3.181617490309819114e-05
56+
-1.021455640559948419e+00
57+
0.000000000000000000e+00
58+
3.551384812675090132e-05
59+
-1.021633100570899177e+00
60+
0.000000000000000000e+00
61+
7.657605381684069333e-05
62+
-1.021837513095323446e+00
63+
0.000000000000000000e+00
64+
5.741698845149804085e-05
65+
-1.022083073542861653e+00
66+
0.000000000000000000e+00
67+
3.051019689016866974e-05
68+
-1.022284054406940923e+00
69+
0.000000000000000000e+00
70+
1.160920758029189404e-05
71+
-1.022354263954097009e+00
72+
0.000000000000000000e+00
73+
6.853375242250866811e-04
74+
-1.020726553087415001e+00
75+
0.000000000000000000e+00
76+
-2.476750464487790571e-04
77+
-1.021164846055963471e+00
78+
0.000000000000000000e+00
79+
1.042038725893374611e-04
80+
-1.021477845759412473e+00
81+
0.000000000000000000e+00
82+
1.078132612144116386e-04
83+
-1.021633493270543536e+00
84+
0.000000000000000000e+00
85+
1.898047742370184514e-04
86+
-1.021832338471513335e+00
87+
0.000000000000000000e+00
88+
1.605551358223333846e-04
89+
-1.022031086967708280e+00
90+
0.000000000000000000e+00
91+
9.463992637281622503e-05
92+
-1.022234850131815254e+00
93+
0.000000000000000000e+00
94+
3.847418450690985012e-05
95+
-1.022309239674765724e+00
96+
0.000000000000000000e+00

0 commit comments

Comments
 (0)