From d9212ed076e0d5c8803b781eb0fda2a5be2a5096 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Tue, 28 Nov 2023 14:03:46 -0600 Subject: [PATCH 01/18] Modifications to try to train KM models --- pyoptmat/chunktime.py | 35 ++++++++++++++++++++++++++++++++--- pyoptmat/models.py | 2 +- pyoptmat/ode.py | 2 +- 3 files changed, 34 insertions(+), 5 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index a84e0c1..bd1c161 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -48,9 +48,8 @@ def newton_raphson_chunk( i = 0 while (i < miter) and torch.any(nR > atol) and torch.any(nR / nR0 > rtol): - x -= solver.solve(J, R) - R, J = fn(x) - nR = torch.norm(R, dim=-1) + dx = solver.solve(J, R) + x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R) i += 1 if i == miter: @@ -60,6 +59,36 @@ def newton_raphson_chunk( return x +def chunk_linesearch(x, dx, fn, R0, sigma = 2.0, c=1e-3, miter = 10): + """ + Backtracking linesearch for the chunk NR algorithm. + + Terminates when the Armijo criteria is reached, or you exceed some maximum iterations + + Args: + x (torch.tensor): initial point + dx (torch.tensor): direction + R0 (torch.tensor): initial residual + + Keyword Args: + sigma (scalar): decrease factor, i.e. alpha /= sigma + c (scalar): stopping criteria + miter (scalar): maximum iterations + """ + alpha = 1.0 + nR0 = torch.norm(R0, dim = -1) + i = 0 + while True: + R, J = fn(x - dx * alpha) + nR = torch.max(torch.norm(R, dim = -1)**2.0) + crit = torch.max(nR0**2.0 + 2.0 * c * alpha * torch.einsum('...i,...i', R0, dx)) + i += 1 + if nR <= crit or i >= miter: + break + alpha /= sigma + return x - dx * alpha, R, J, torch.max(torch.norm(R, dim = -1)), alpha + + class BidiagonalOperator(torch.nn.Module): """ diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 19d4d1d..a8d9e6a 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -545,7 +545,7 @@ def RJ(erate): return R[..., None], J[..., None, None] - erate, _ = solvers.newton_raphson(RJ, erate_guess) + erate, _ = solvers.newton_raphson(RJ, erate_guess, atol = 1.0e-2) yp = y.clone() yp[..., 0] = cs ydot, J, Je, _ = self.model(t, yp, erate[..., 0], cT) diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index cb7989f..0b0b946 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -293,7 +293,7 @@ def __init__( block_size=1, rtol=1.0e-6, atol=1.0e-8, - miter=100, + miter=200, linear_solve_method="direct", direct_solve_method="thomas", direct_solve_min_size=0, From 03b7d63d1c0b1ec3fa4b347dd85393f16a97eaf0 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 30 Nov 2023 12:49:22 -0600 Subject: [PATCH 02/18] Continued tweaks to integration routines --- pyoptmat/chunktime.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index bd1c161..70e9629 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -47,7 +47,7 @@ def newton_raphson_chunk( nR0 = nR i = 0 - while (i < miter) and torch.any(nR > atol) and torch.any(nR / nR0 > rtol): + while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))): dx = solver.solve(J, R) x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R) i += 1 @@ -86,7 +86,7 @@ def chunk_linesearch(x, dx, fn, R0, sigma = 2.0, c=1e-3, miter = 10): if nR <= crit or i >= miter: break alpha /= sigma - return x - dx * alpha, R, J, torch.max(torch.norm(R, dim = -1)), alpha + return x - dx * alpha, R, J, torch.norm(R, dim = -1), alpha From d00999e990b062b3eec8a1060c01651545c49a1e Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 30 Nov 2023 14:16:16 -0600 Subject: [PATCH 03/18] Added generalized chunk sizes --- pyoptmat/chunktime.py | 12 ++++++-- pyoptmat/ode.py | 69 ++++++++++++++++++++++++++++++------------- test/test_ode.py | 43 +++++++++++++++++++++++++++ 3 files changed, 102 insertions(+), 22 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index 70e9629..7187cb5 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -20,7 +20,8 @@ def newton_raphson_chunk( - fn, x0, solver, rtol=1e-6, atol=1e-10, miter=100, throw_on_fail=False + fn, x0, solver, rtol=1e-6, atol=1e-10, miter=100, throw_on_fail=False, + linesearch = False ): """ Solve a nonlinear system with Newton's method with a tensor for a @@ -36,6 +37,7 @@ def newton_raphson_chunk( atol (float): nonlinear absolute tolerance miter (int): maximum number of nonlinear iterations throw_on_fail (bool): throw exception if solve fails + linesearch (bool): if true use backtracking linesearch Returns: torch.tensor: solution to system of equations @@ -49,7 +51,13 @@ def newton_raphson_chunk( while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))): dx = solver.solve(J, R) - x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R) + if linesearch: + x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R) + else: + x -= dx + R, J = fn(x) + nR = torch.norm(R, dim = -1) + print(i, torch.max(nR)) i += 1 if i == miter: diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index 0b0b946..5946489 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -275,6 +275,7 @@ class FixedGridBlockSolver: rtol (float): relative tolerance for Newton's method atol (float): absolute tolerance for Newton's method miter (int): maximum number of Newton iterations + linesearch (bool): whether to use linear_solve_method (str): method to solve batched sparse Ax = b, options are currently "direct" or "dense" direct_solve_method (str): method to use for the direct solver, options are @@ -282,7 +283,8 @@ class FixedGridBlockSolver: direct_solve_min_size (int): minimum PCR block size for the hybrid approach adjoint_params: parameters to track for the adjoint backward pass guess_type (string): strategy for initial guess, options are "zero" and "previous" - throw_on_fail (bool): if true throw an exception if the implicit solve fails + throw_on_fail (bool): if true throw an exception if the implicit solve fails, + offset_step (int): use a special, smaller chunk size for the first step """ def __init__( @@ -294,12 +296,14 @@ def __init__( rtol=1.0e-6, atol=1.0e-8, miter=200, + linesearch = False, linear_solve_method="direct", direct_solve_method="thomas", direct_solve_min_size=0, adjoint_params=None, guess_type="zero", throw_on_fail=False, + offset_step = 0, **kwargs, ): # Store basic info about the system @@ -318,6 +322,7 @@ def __init__( self.rtol = rtol self.atol = atol self.miter = miter + self.linesearch = linesearch # Store for later self.adjoint_params = adjoint_params @@ -347,6 +352,9 @@ def __init__( # Throw exception on failed solve self.throw_on_fail = throw_on_fail + # Special first chunk size + self.offset_step = offset_step + # Cached solutions self.t = None self.result = None @@ -370,14 +378,15 @@ def integrate(self, t, cache_adjoint=False): t.shape[0], *self.y0.shape, dtype=self.y0.dtype, device=self.y0.device ) result[0] = self.y0 + incs = self._gen_increments(t) - for k in range(1, t.shape[0], self.n): - result[k : k + self.n] = self.block_update( - t[k : k + self.n], - t[k - 1], - result[k - 1], + for k1,k2 in zip(incs[:-1], incs[1:]): + result[k1 : k2] = self.block_update( + t[k1 : k2], + t[k1 - 1], + result[k1 - 1], self.func, - self._initial_guess(result, k), + self._initial_guess(result, k1, k2-k1), ) # Store for the backward pass, if we're going to do that @@ -387,22 +396,38 @@ def integrate(self, t, cache_adjoint=False): return result - def _initial_guess(self, result, k): + def _gen_increments(self, t): + """ + Generate the increments in time to use the chunk integrate the equations + + Args: + t (torch.tensor): timesteps requested + """ + steps = [1] + ntotal = t.shape[0] + if self.offset_step > 0: + steps += [self.offset_step + 1] + steps += list(range(steps[-1],ntotal, self.n))[1:] + [ntotal] + + return steps + + def _initial_guess(self, result, k, nchunk): """ Form the initial guess Args: result (torch.tensor): currently-populated results k (int): current time step + nchunk (int): current chunk size """ if self.guess_type == "zero": - guess = torch.zeros_like(result[k : k + self.n]) + guess = torch.zeros_like(result[k : k + nchunk]) elif self.guess_type == "previous": - if k - self.n - 1 < 0: - guess = torch.zeros_like(result[k : k + self.n]) + if k - nchunk - 1 < 0: + guess = torch.zeros_like(result[k : k + nchunk]) else: - guess = result[(k - self.n) : k] - result[k - self.n - 1].unsqueeze(0) - blk = self.n - result[k : k + self.n].shape[0] + guess = result[(k - nchunk) : k] - result[k - nchunk - 1].unsqueeze(0) + blk = nchunk - result[k : k + nchunk].shape[0] guess = guess[blk:] else: raise ValueError(f"Unknown initial guess strategy {self.guess_type}!") @@ -439,18 +464,21 @@ def rewind(self, output_grad): # Calculate starts at last gradient prev_adjoint = output_grad[0] + + # Generate reverse increments + incs = [1 + self.t.shape[0] - r for r in self._gen_increments(self.t)[::-1]] - for k in range(1, self.t.shape[0], self.n): + for k1, k2 in zip(incs[:-1], incs[1:]): # Could also cache these of course _, J = self.func( - self.t[k - 1 : k + self.n], self.result[k - 1 : k + self.n] + self.t[k1 - 1 : k2], self.result[k1 - 1 : k2] ) full_adjoint = self.scheme.update_adjoint( - self.t[k - 1 : k + self.n].diff(dim=0), + self.t[k1 - 1 : k2].diff(dim=0), J, prev_adjoint, - output_grad[k - 1 : k + self.n], + output_grad[k1 - 1 : k2], solver=self.direct_solver, ) @@ -460,10 +488,10 @@ def rewind(self, output_grad): if len(self.adjoint_params) > 0: grad_result = self.scheme.accumulate( grad_result, - self.t[k - 1 : k + self.n], - self.result[k - 1 : k + self.n], + self.t[k1 - 1 : k2], + self.result[k1 - 1 : k2], full_adjoint, - output_grad[k - 1 : k + self.n], + output_grad[k1 - 1 : k2], self._get_param_partial, ) @@ -512,6 +540,7 @@ def RJ(dy): atol=self.atol, miter=self.miter, throw_on_fail=self.throw_on_fail, + linesearch = self.linesearch ) return dy.reshape(self.batch_size, n, self.prob_size).transpose( diff --git a/test/test_ode.py b/test/test_ode.py index e433921..399ebb1 100644 --- a/test/test_ode.py +++ b/test/test_ode.py @@ -181,6 +181,49 @@ def test_backward(self): self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) +class TestFallingBatchUnequal(unittest.TestCase): + def setUp(self): + self.m = 1.0 / 4.0 + self.w = 8.0 + self.k = 2.0 + + self.nbatch = 20 + self.y0 = torch.zeros(self.nbatch, 2) + + self.nsteps = 1000 + self.nchunk = 10 + self.times = torch.tensor( + np.array([np.linspace(0, 5.0, self.nsteps) for i in range(self.nbatch)]).T + ) + + self.model = FallingODE(self.m, self.w, self.k) + + def test_forward(self): + exact = self.model.exact(self.times, self.y0) + numerical = ode.odeint( + self.model, + self.y0, + self.times, + method="forward-euler", + block_size=self.nchunk, + offset_step = 2 + ) + + self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) + + def test_backward(self): + exact = self.model.exact(self.times, self.y0) + numerical = ode.odeint( + self.model, + self.y0, + self.times, + method="backward-euler", + block_size=self.nchunk, + offset_step = 2 + ) + + self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) + class FallingParameterizedODE(torch.nn.Module): """ From 205ff6ac2ecc5ab03ab2ad4fe73b63397ffc6217 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 30 Nov 2023 17:10:42 -0600 Subject: [PATCH 04/18] More progress --- .../creep-fatigue/statistical/infer.py | 1 + .../tension/statistical/infer.py | 2 +- pyoptmat/chunktime.py | 1 - pyoptmat/models.py | 33 +++++--- pyoptmat/ode.py | 8 +- pyoptmat/optimize.py | 47 +++++++++-- pyoptmat/solvers.py | 79 +++++++++++++++++++ 7 files changed, 150 insertions(+), 21 deletions(-) diff --git a/examples/structural-inference/creep-fatigue/statistical/infer.py b/examples/structural-inference/creep-fatigue/statistical/infer.py index 2359ce7..a7f28dd 100755 --- a/examples/structural-inference/creep-fatigue/statistical/infer.py +++ b/examples/structural-inference/creep-fatigue/statistical/infer.py @@ -94,6 +94,7 @@ def make(n, eta, s0, R, d, C, g, **kwargs): scale_scale_priors, eps, include_noise=False, + use_cached_guess= True ).to(device) # 4) Get the guide diff --git a/examples/structural-inference/tension/statistical/infer.py b/examples/structural-inference/tension/statistical/infer.py index 9a015d3..a998c41 100755 --- a/examples/structural-inference/tension/statistical/infer.py +++ b/examples/structural-inference/tension/statistical/infer.py @@ -91,7 +91,7 @@ def make(n, eta, s0, R, d, **kwargs): model = optimize.HierarchicalStatisticalModel( lambda *args, **kwargs: make(*args, block_size = time_chunk_size, direct_solve_method = linear_solve_method, **kwargs), - names, loc_loc_priors, loc_scale_priors, scale_scale_priors, eps + names, loc_loc_priors, loc_scale_priors, scale_scale_priors, eps, use_cached_guess = False ).to(device) # 4) Get the guide diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index 7187cb5..f42beb5 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -57,7 +57,6 @@ def newton_raphson_chunk( x -= dx R, J = fn(x) nR = torch.norm(R, dim = -1) - print(i, torch.max(nR)) i += 1 if i == miter: diff --git a/pyoptmat/models.py b/pyoptmat/models.py index a8d9e6a..3cbedc1 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -283,7 +283,7 @@ class ModelIntegrator(nn.Module): """ - def __init__(self, model, *args, use_adjoint=True, **kwargs): + def __init__(self, model, *args, use_adjoint=True, bisect_first = False, **kwargs): super().__init__(*args) self.model = model self.use_adjoint = use_adjoint @@ -294,6 +294,8 @@ def __init__(self, model, *args, use_adjoint=True, **kwargs): else: self.imethod = ode.odeint + self.bisect_first = bisect_first + def solve_both(self, times, temperatures, idata, control): """ Solve for either strain or stress control at once @@ -327,6 +329,7 @@ def solve_both(self, times, temperatures, idata, control): base_interpolator, temperature_interpolator, control, + bisect_first = self.bisect_first ) return self.imethod(bmodel, init, times, **self.kwargs_for_integration) @@ -435,7 +438,7 @@ class BothBasedModel(nn.Module): indices: split into strain and stress control """ - def __init__(self, model, rate_fn, base_fn, T_fn, control, *args, **kwargs): + def __init__(self, model, rate_fn, base_fn, T_fn, control, bisect_first = False, *args, **kwargs): super().__init__(*args, **kwargs) self.model = model self.rate_fn = rate_fn @@ -445,7 +448,7 @@ def __init__(self, model, rate_fn, base_fn, T_fn, control, *args, **kwargs): self.emodel = StrainBasedModel(self.model, self.rate_fn, self.T_fn) self.smodel = StressBasedModel( - self.model, self.rate_fn, self.base_fn, self.T_fn + self.model, self.rate_fn, self.base_fn, self.T_fn, bisect_first = bisect_first ) def forward(self, t, y): @@ -514,12 +517,15 @@ class StressBasedModel(nn.Module): T_fn: T(t) """ - def __init__(self, model, srate_fn, stress_fn, T_fn, *args, **kwargs): + def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate = 1e3, bisect_first = False, *args, **kwargs): super().__init__(*args, **kwargs) self.model = model self.srate_fn = srate_fn self.stress_fn = stress_fn self.T_fn = T_fn + self.min_erate = min_erate + self.max_erate = max_erate + self.bisect_first = bisect_first def forward(self, t, y): """ @@ -533,28 +539,33 @@ def forward(self, t, y): cs = self.stress_fn(t) cT = self.T_fn(t) - erate_guess = torch.zeros_like(y[..., 0])[..., None] - def RJ(erate): yp = y.clone() yp[..., 0] = cs - ydot, _, Je, _ = self.model(t, yp, erate[..., 0], cT) + ydot, _, Je, _ = self.model(t, yp, erate, cT) R = ydot[..., 0] - csr J = Je[..., 0] - return R[..., None], J[..., None, None] + return R, J + + if self.bisect_first: + erate = solvers.scalar_bisection_newton(RJ, + torch.ones_like(y[...,0]) * self.min_erate, + torch.ones_like(y[...,0]) * self.max_erate) + else: + erate = solvers.scalar_newton(RJ, + torch.zeros_like(y[...,0])) - erate, _ = solvers.newton_raphson(RJ, erate_guess, atol = 1.0e-2) yp = y.clone() yp[..., 0] = cs - ydot, J, Je, _ = self.model(t, yp, erate[..., 0], cT) + ydot, J, Je, _ = self.model(t, yp, erate, cT) # Rescale the jacobian J[..., 0, :] = -J[..., 0, :] / Je[..., 0][..., None] J[..., :, 0] = 0 # Insert the strain rate - ydot[..., 0] = erate[..., 0] + ydot[..., 0] = erate return ydot, J diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index 5946489..007de21 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -283,6 +283,7 @@ class FixedGridBlockSolver: direct_solve_min_size (int): minimum PCR block size for the hybrid approach adjoint_params: parameters to track for the adjoint backward pass guess_type (string): strategy for initial guess, options are "zero" and "previous" + guess_history (torch.tensor): complete load history used for guess, overrides guess_type throw_on_fail (bool): if true throw an exception if the implicit solve fails, offset_step (int): use a special, smaller chunk size for the first step """ @@ -302,6 +303,7 @@ def __init__( direct_solve_min_size=0, adjoint_params=None, guess_type="zero", + guess_history=None, throw_on_fail=False, offset_step = 0, **kwargs, @@ -348,6 +350,7 @@ def __init__( # Initial guess for integration self.guess_type = guess_type + self.guess_history = guess_history # Throw exception on failed solve self.throw_on_fail = throw_on_fail @@ -380,6 +383,7 @@ def integrate(self, t, cache_adjoint=False): result[0] = self.y0 incs = self._gen_increments(t) + for k1,k2 in zip(incs[:-1], incs[1:]): result[k1 : k2] = self.block_update( t[k1 : k2], @@ -420,7 +424,9 @@ def _initial_guess(self, result, k, nchunk): k (int): current time step nchunk (int): current chunk size """ - if self.guess_type == "zero": + if self.guess_history is not None: + guess = self.guess_history[k : k + nchunk] - result[k - 1] + elif self.guess_type == "zero": guess = torch.zeros_like(result[k : k + nchunk]) elif self.guess_type == "previous": if k - nchunk - 1 < 0: diff --git a/pyoptmat/optimize.py b/pyoptmat/optimize.py index c9d80e9..26a88bc 100644 --- a/pyoptmat/optimize.py +++ b/pyoptmat/optimize.py @@ -124,7 +124,7 @@ class DeterministicModel(Module): ics (list(torch.tensor)): initial conditions to use for each parameter """ - def __init__(self, maker, names, ics): + def __init__(self, maker, names, ics, use_cached_guess = False): super().__init__() self.maker = maker @@ -136,6 +136,9 @@ def __init__(self, maker, names, ics): for name, ic in zip(names, ics): setattr(self, name, Parameter(ic)) + self.use_cached_guess = use_cached_guess + self.cached_solution = None + def get_params(self): """ Return the parameters for input to the model @@ -155,12 +158,18 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control): exp_types (torch.tensor): experiment types, as integers exp_control (torch.tensor): stress/strain control flag """ - model = self.maker(*self.get_params()) + if self.use_cached_guess: + model = self.maker(*self.get_params(), guess_history = self.cached_solution) + else: + model = self.maker(*self.get_params()) predictions = model.solve_both( exp_data[0], exp_data[1], exp_data[2], exp_control ) + if self.use_cached_guess: + self.cached_solution = predictions.detach().clone() + return experiments.convert_results(predictions[:, :, 0], exp_cycles, exp_types) @@ -189,7 +198,7 @@ class StatisticalModel(PyroModule): entry i represents the noise in test type i """ - def __init__(self, maker, names, locs, scales, eps, nan_num=False): + def __init__(self, maker, names, locs, scales, eps, nan_num=False, use_cached_guess = False): super().__init__() self.maker = maker @@ -205,6 +214,9 @@ def __init__(self, maker, names, locs, scales, eps, nan_num=False): self.nan_num = nan_num + self.use_cached_guess = use_cached_guess + self.cached_solution = None + def get_params(self): """ Return the sampled parameters for input to the model @@ -229,10 +241,17 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control, exp_results=None Keyword Args: exp_results (torch.tensor): true results for conditioning """ - model = self.maker(*self.get_params()) + if self.use_cached_guess: + model = self.maker(*self.get_params(), guess_history = self.cached_solution) + else: + model = self.maker(*self.get_params()) + predictions = model.solve_both( exp_data[0], exp_data[1], exp_data[2], exp_control ) + if self.use_cached_guess: + self.cached_solution = predictions.detach().clone() + results = experiments.convert_results( predictions[:, :, 0], exp_cycles, exp_types ) @@ -318,6 +337,7 @@ def __init__( param_suffix="_param", include_noise=False, weights=None, + use_cached_guess = False ): super().__init__() @@ -395,6 +415,9 @@ def __init__( # This annoyance is required to make the adjoint solver work self.extra_param_names = [] + self.use_cached_guess = use_cached_guess + self.cached_solution = None + @property def nparams(self): """ @@ -531,13 +554,23 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control, exp_results=None scale=self._make_weight_tensor(exp_types) ): # Sample the bottom level parameters - bmodel = self.maker( - *self.sample_bot(), extra_params=self.get_extra_params() - ) + if self.use_cached_guess: + bmodel = self.maker( + *self.sample_bot(), extra_params=self.get_extra_params(), + guess_history = self.cached_solution + ) + else: + bmodel = self.maker( + *self.sample_bot(), extra_params=self.get_extra_params() + ) + # Generate the results predictions = bmodel.solve_both( exp_data[0], exp_data[1], exp_data[2], exp_control ) + if self.use_cached_guess: + self.cached_solution = predictions.detach().clone() + # Process the results results = experiments.convert_results( predictions[:, :, 0], exp_cycles, exp_types diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index 5044715..a6014c2 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -7,6 +7,85 @@ import warnings import torch +def scalar_bisection(fn, a, b, atol = 1.0e-6, miter = 100): + """ + Solve logically scalar equations with bisection + + Args: + fn (function): function returning scalar residual and jacobian + a (torch.tensor): lower bound + b (torch.tensor): upper bound + + Keyword Args: + atol : absolute tolerance for convergence + miter (int): max number of iterations + """ + Ra, _ = fn(a) + Rb, _ = fn(b) + + if not torch.all((torch.sign(Ra) + torch.sign(Rb)) == 0): + raise RuntimeError("Initial values do not bisect in bisection solver") + + c = (a+b) / 2.0 + Rc, _ = fn(c) + + for i in range(miter): + if torch.all(torch.abs(Rc) < atol): + break + + ac = torch.sign(Ra) == torch.sign(Rc) + bc = torch.sign(Rb) == torch.sign(Rc) + a[ac] = c[ac] + b[bc] = c[bc] + + c = (a+b) / 2.0 + Rc, _ = fn((a+b) / 2.0) + + return c + +def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100): + """ + Solve logically scalar equations with Newton's method + + Args: + fn (function): function returning scalar residual and jacobian + x0 (torch.tensor): initial guess + + Keyword Args: + atol (float): absolute tolerance for convergence + miter (int): maximum number of iterations + """ + x = x0 + R, J = fn(x) + + for i in range(miter): + if torch.all(torch.abs(R) < atol): + break + + x -= R / J + + R, J = fn(x) + else: + warnings.warn("Scalar implicit solve did not succeed. Results may be inaccurate...") + + return x + +def scalar_bisection_newton(fn, a, b, atol = 1.0e-6, miter = 100, biter = 10): + """ + Solve logically scalar equations by switching from bisection to Newton's method + + Args: + fn (function): function returning scalar residual and jacobian + a (torch.tensor): lower bound + b (torch.tensor): upper bound + + Keyword Args: + atol : absolute tolerance for convergence + biter: initial number of bisection iterations + miter (int): max number of iterations for Newton's method + """ + x = scalar_bisection(fn, a, b, atol = atol, miter = biter) + return scalar_newton(fn, x, atol = atol, miter = miter) def newton_raphson_bt( fn, x0, linsolver="lu", rtol=1e-6, atol=1e-10, miter=100, max_bt=5 From 46d05f39101f65cad4937341c40b488442b9ba39 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Tue, 2 Jan 2024 16:22:42 -0600 Subject: [PATCH 05/18] Fixed in place operation error --- pyoptmat/models.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 3cbedc1..8625289 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -540,15 +540,16 @@ def forward(self, t, y): cT = self.T_fn(t) def RJ(erate): - yp = y.clone() - yp[..., 0] = cs - ydot, _, Je, _ = self.model(t, yp, erate, cT) + ydot, _, Je, _ = self.model(t, + torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), + erate.detach(), cT) R = ydot[..., 0] - csr J = Je[..., 0] return R, J + # Doing the detach here actually makes the parameter gradients wrong... if self.bisect_first: erate = solvers.scalar_bisection_newton(RJ, torch.ones_like(y[...,0]) * self.min_erate, @@ -557,15 +558,13 @@ def RJ(erate): erate = solvers.scalar_newton(RJ, torch.zeros_like(y[...,0])) - yp = y.clone() - yp[..., 0] = cs - ydot, J, Je, _ = self.model(t, yp, erate, cT) + ydot, J, Je, _ = self.model(t, + torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), + erate, cT) - # Rescale the jacobian - J[..., 0, :] = -J[..., 0, :] / Je[..., 0][..., None] - J[..., :, 0] = 0 + # Corrected jacobian + row1 = torch.cat([torch.zeros_like(J[...,0,0]).unsqueeze(-1), -J[..., 0, 1:] / Je[...,0][...,None]], dim = -1).unsqueeze(-2) + rest = torch.cat([torch.zeros_like(J[...,1:,0]).unsqueeze(-1), J[...,1:,1:]], dim = -1) + jac = torch.cat([row1,rest], dim = -2) - # Insert the strain rate - ydot[..., 0] = erate - - return ydot, J + return torch.cat([erate.unsqueeze(-1), ydot[...,1:]], dim = -1), jac From ffa368668c18ca31b08c3d6d176deef68bc6f129 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 4 Jan 2024 08:35:00 -0600 Subject: [PATCH 06/18] Tweaks to speed things up --- pyoptmat/chunktime.py | 23 ++++++++++-------- pyoptmat/models.py | 54 ++++++++++++++++++++++--------------------- 2 files changed, 41 insertions(+), 36 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index f42beb5..b7bdb6f 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -52,7 +52,7 @@ def newton_raphson_chunk( while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))): dx = solver.solve(J, R) if linesearch: - x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R) + x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R, rtol, atol) else: x -= dx R, J = fn(x) @@ -66,34 +66,37 @@ def newton_raphson_chunk( return x -def chunk_linesearch(x, dx, fn, R0, sigma = 2.0, c=1e-3, miter = 10): +def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=1e-3, miter = 10): """ Backtracking linesearch for the chunk NR algorithm. - Terminates when the Armijo criteria is reached, or you exceed some maximum iterations + Terminates when the Armijo criteria is reached, or you exceed some maximum iterations, or when you would meet the + convergence requirements with the current alpha Args: x (torch.tensor): initial point dx (torch.tensor): direction R0 (torch.tensor): initial residual + overall_rtol (scalar): Newton relative tolerance + overall_atol (scalar): Newton absolute tolerance Keyword Args: sigma (scalar): decrease factor, i.e. alpha /= sigma c (scalar): stopping criteria miter (scalar): maximum iterations """ - alpha = 1.0 + alpha = torch.ones(x.shape[:-1], device = x.device) nR0 = torch.norm(R0, dim = -1) i = 0 while True: - R, J = fn(x - dx * alpha) - nR = torch.max(torch.norm(R, dim = -1)**2.0) - crit = torch.max(nR0**2.0 + 2.0 * c * alpha * torch.einsum('...i,...i', R0, dx)) + R, J = fn(x - dx * alpha.unsqueeze(-1)) + nR = torch.norm(R, dim = -1)**2.0 + crit = nR0**2.0 + 2.0 * c * alpha * torch.einsum('...i,...i', R0, dx) i += 1 - if nR <= crit or i >= miter: + if torch.all(nR <= crit) or i >= miter or torch.all(torch.sqrt(nR) < overall_atol) or torch.all(torch.sqrt(nR)/nR0 < overall_rtol): break - alpha /= sigma - return x - dx * alpha, R, J, torch.norm(R, dim = -1), alpha + alpha[nR > crit] /= sigma + return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim = -1), alpha diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 8625289..a56054e 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -315,19 +315,14 @@ def solve_both(self, times, temperatures, idata, control): # Likely if this happens dt = 0 rates[torch.isnan(rates)] = 0 - rate_interpolator = utility.ArbitraryBatchTimeSeriesInterpolator(times, rates) - base_interpolator = utility.ArbitraryBatchTimeSeriesInterpolator(times, idata) - temperature_interpolator = utility.ArbitraryBatchTimeSeriesInterpolator( - times, temperatures - ) - init = torch.zeros(times.shape[1], self.model.nhist, device=idata.device) bmodel = BothBasedModel( self.model, - rate_interpolator, - base_interpolator, - temperature_interpolator, + times, + rates, + idata, + temperatures, control, bisect_first = self.bisect_first ) @@ -438,17 +433,23 @@ class BothBasedModel(nn.Module): indices: split into strain and stress control """ - def __init__(self, model, rate_fn, base_fn, T_fn, control, bisect_first = False, *args, **kwargs): + def __init__(self, model, times, rates, base, temps, control, bisect_first = False, *args, **kwargs): super().__init__(*args, **kwargs) self.model = model - self.rate_fn = rate_fn - self.base_fn = base_fn - self.T_fn = T_fn self.control = control - self.emodel = StrainBasedModel(self.model, self.rate_fn, self.T_fn) + self.econtrol = self.control == 0 + self.scontrol = self.control == 1 + + self.emodel = StrainBasedModel(self.model, + utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.econtrol], rates[...,self.econtrol]), + utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.econtrol], temps[...,self.econtrol])) self.smodel = StressBasedModel( - self.model, self.rate_fn, self.base_fn, self.T_fn, bisect_first = bisect_first + self.model, + utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], rates[...,self.scontrol]), + utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], base[...,self.scontrol]), + utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], temps[...,self.scontrol]), + bisect_first = bisect_first ) def forward(self, t, y): @@ -460,20 +461,21 @@ def forward(self, t, y): t: input times y: input state """ - strain_rates, strain_jacs = self.emodel(t, y) - stress_rates, stress_jacs = self.smodel(t, y) - - actual_rates = torch.zeros_like(strain_rates) + n = (y.shape[-1],) + base = y.shape[:-1] - e_control = self.control == 0 - s_control = self.control == 1 + actual_rates = torch.zeros(base + n, device = t.device) + actual_jacs = torch.zeros(base + n + n, device = t.device) - actual_rates[..., e_control, :] = strain_rates[..., e_control, :] - actual_rates[..., s_control, :] = stress_rates[..., s_control, :] + if torch.any(self.econtrol): + strain_rates, strain_jacs = self.emodel(t[...,self.econtrol], y[...,self.econtrol,:]) + actual_rates[...,self.econtrol,:] = strain_rates + actual_jacs[...,self.econtrol,:,:] = strain_jacs - actual_jacs = torch.zeros_like(strain_jacs) - actual_jacs[..., e_control, :, :] = strain_jacs[..., e_control, :, :] - actual_jacs[..., s_control, :, :] = stress_jacs[..., s_control, :, :] + if torch.any(self.scontrol): + stress_rates, stress_jacs = self.smodel(t[...,self.scontrol], y[...,self.scontrol,:]) + actual_rates[...,self.scontrol,:] = stress_rates + actual_jacs[...,self.scontrol,:,:] = stress_jacs return actual_rates, actual_jacs From 76f28136f32f9a11ed1435f6390f22545fd67db2 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Fri, 5 Jan 2024 17:47:14 -0600 Subject: [PATCH 07/18] Sigh, still issues --- pyoptmat/flowrules.py | 419 +++++++++++++++++++++++++++++++++++++++++ test/test_flowrules.py | 130 +++++++++++-- 2 files changed, 537 insertions(+), 12 deletions(-) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index 249ba36..a3b2dd3 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -1125,6 +1125,95 @@ def history_rate(self, s, h, t, T, e): """ return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) +class ApproximatePerfectRateIndependentPlasticity(FlowRule): + """ + Perfect viscoplasticity defined as + + .. math:: + + \\dot{\\varepsilon}_{in}=\\left(\\frac{\\left|\\sigma\\right|}{\\eta}\\right)^{n}\\operatorname{sign}\\left(\\sigma\\right) + + \\dot{h} = \\emptyset + + Args: + n (|TP|): rate sensitivity + eta (|TP|): flow viscosity + """ + + def __init__(self, sy, n = 20.0): + super().__init__() + self.sy = sy + self.n = n + + def flow_rate(self, s, h, t, T, e): + """ + The uniaxial flow rate itself and the derivative + with respect to stress + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + tuple(torch.tensor, torch.tensor): the flow rate and the derivative + of the flow rate with + respect to stress + """ + return ( + torch.abs(e) * utility.macaulay(torch.abs(s) - self.sy(T)) ** self.n * torch.sign(s), + torch.abs(e) * self.n * utility.macaulay(torch.abs(s) - self.sy(T))**(self.n-1) * utility.heaviside(torch.abs(s) - self.sy(T)) * torch.sign(s) + ) + + def dflow_derate(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the total strain rate + + The superclass implementation provides a default of zero with the + right shape. + + Args: + s (torch.tensor): stress + h (torch.tensor): internal variables + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: derivative of flow rate with respect to the + internal variables + """ + return torch.sign(e) * utility.macaulay(torch.abs(s) - self.sy(T)) ** self.n * torch.sign(s) + + @property + def nhist(self): + """ + The number of internal variables + + Here 0... + """ + return 0 + + def history_rate(self, s, h, t, T, e): + """ + The history rate and the derivative of the history rate with respect + to the current history + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + + Returns: + tuple(torch.tensor, torch.tensor): the history rate and the + derivative of the history rate + with respect to history + """ + return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) + class IsoKinViscoplasticity(FlowRule): """ @@ -1425,3 +1514,333 @@ def dhist_derate(self, s, h, t, T, e): ], dim=-1, ) + + +class ApproximateIsoKinRateIndependentPlasticity(FlowRule): + """ + Viscoplasticity with isotropic and kinematic hardening, defined as + + .. math:: + + \\dot{\\varepsilon}_{in}=\\left\\langle \\frac{\\left|\\sigma-x\\right|-s_{0}-k}{\\eta}\\right\\rangle ^{n}\\operatorname{sign}\\left(\\sigma-X\\right) + + and where the :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and + :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects determine the + history rate. + + The :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and + :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects each define both a + set of internal variables, including the corresponding rate forms and Jacobians, + but also a map from those internal variables to the isotropic hardening + value :math:`k` (for :py:class:`pyoptmat.hardening.IsotropicHardeningModel`) + and the kinematic hardening value :math:`x` + (for :py:class:`pyoptmat.hardening.KinematicHardeningModel`), along with + the derivatives of those maps. All this information is required to assemble + the information this class needs to provide. + + Args: + n (|TP|): rate sensitivity + eta (|TP|): flow viscosity + s0 (|TP|): initial value of flow stress (i.e. the threshold stress) + isotropic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the isotropic hardening model + kinematic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the kinematic hardening model + """ + + def __init__(self, sy, isotropic, kinematic, n = 20.0): + super().__init__() + self.isotropic = isotropic + self.kinematic = kinematic + self.n = n + self.sy = sy + + def flow_rate(self, s, h, t, T, e): + """ + The flow rate itself and the derivative with respect to stress + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + tuple(torch.tensor, torch.tensor): the flow rate and the derivative of the flow rate with + respect to stress + """ + ih = self.isotropic.value(h[..., : self.isotropic.nhist]) + kh = self.kinematic.value(h[..., self.isotropic.nhist :]) + + return ( + torch.abs(e) * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) + ** self.n + * torch.sign(s - kh), + torch.abs(e) * self.n + * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) + ** (self.n - 1), + ) + + def dflow_diso(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the isotropic hardening + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: the derivative of the flow rate + """ + ih = self.isotropic.value(h[..., : self.isotropic.nhist]) + kh = self.kinematic.value( + h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] + ) + + iv = (torch.abs(s - kh) - self.sy(T) - ih) + + return ( + -torch.abs(e) * self.n + * utility.macaulay(iv) ** (self.n - 1) + * torch.sign(s - kh) + ) + + def dflow_dkin(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the kinematic hardening + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: the derivative of the flow rate + """ + ih = self.isotropic.value(h[..., : self.isotropic.nhist]) + kh = self.kinematic.value( + h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] + ) + + return ( + -torch.abs(e) * self.n + * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) + ** (self.n - 1) + ) + + def dflow_derate(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the total strain rate + + The superclass implementation provides a default of zero with the + right shape. + + Args: + s (torch.tensor): stress + h (torch.tensor): internal variables + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: derivative of flow rate with respect to the + internal variables + """ + ih = self.isotropic.value(h[..., : self.isotropic.nhist]) + kh = self.kinematic.value( + h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] + ) + + return torch.sign(e) * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih))**self.n * torch.sign(s - kh) + + @property + def nhist(self): + """ + The number of internal variables, here the sum from the isotropic + and kinematic hardening models + """ + return self.isotropic.nhist + self.kinematic.nhist + + def history_rate(self, s, h, t, T, e): + """ + The vector of the rates of the internal variables split into + portions defined by each hardening model + + The first chunk of entries is for the isotropic hardening, + the second for the kinematic hardening. + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + tuple(torch.tensor, torch.tensor): the history rate and the + derivative of the history rate + with respect to history + """ + erate, _ = self.flow_rate(s, h, t, T, e) + + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] + + # History evolution rate + hrate = torch.cat( + [ + self.isotropic.history_rate(s, hiso, t, erate, T, e), + self.kinematic.history_rate(s, hkin, t, erate, T, e), + ], + dim=-1, + ) + + # Jacobian contribution + hdiv = torch.cat( + [ + torch.cat( + [ + self.isotropic.dhistory_rate_dhistory(s, hiso, t, erate, T, e) + + utility.mbmm( + self.isotropic.dhistory_rate_derate( + s, hiso, t, erate, T, e + ), + utility.mbmm( + self.dflow_diso(s, h, t, T, e)[..., None, None], + self.isotropic.dvalue(hiso)[..., None, :], + ), + ), + self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e) + * self.dflow_dkin(s, h, t, T, e)[..., None, None] + * self.kinematic.dvalue(hkin)[..., None, :], + ], + dim=-1, + ), + torch.cat( + [ + utility.mbmm( + self.kinematic.dhistory_rate_derate( + s, hkin, t, erate, T, e + ), + utility.mbmm( + self.dflow_diso(s, h, t, T, e)[..., None, None], + self.isotropic.dvalue(hiso)[..., None, :], + ), + ), + self.kinematic.dhistory_rate_dhistory(s, hkin, t, erate, T, e) + + utility.mbmm( + self.kinematic.dhistory_rate_derate( + s, hkin, t, erate, T, e + ), + utility.mbmm( + self.dflow_dkin(s, h, t, T, e)[..., None, None], + self.kinematic.dvalue(hkin)[..., None, :], + ), + ), + ], + dim=-1, + ), + ], + dim=-2, + ) + + return hrate, hdiv + + def dflow_dhist(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the internal variables + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: the derivative of the flow rate + """ + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] + + return torch.cat( + [ + self.dflow_diso(s, h, t, T, e)[..., None] * self.isotropic.dvalue(hiso), + self.dflow_dkin(s, h, t, T, e)[..., None] * self.kinematic.dvalue(hkin), + ], + dim=-1, + ).unsqueeze(-2) + + def dhist_dstress(self, s, h, t, T, e): + """ + The derivative of the history rate with respect to the stress + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: the derivative of the flow rate + """ + erate, derate = self.flow_rate(s, h, t, T, e) + + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] + + return torch.cat( + [ + self.isotropic.dhistory_rate_dstress(s, hiso, t, erate, T, e) + + utility.mbmm( + self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e), + derate[..., None, None], + )[..., 0], + self.kinematic.dhistory_rate_dstress(s, hkin, t, erate, T, e) + + utility.mbmm( + self.kinematic.dhistory_rate_derate(s, hkin, t, erate, T, e), + derate[..., None, None], + )[..., 0], + ], + dim=-1, + ) + + def dhist_derate(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the total strain rate + + Args: + s (torch.tensor): stress + h (torch.tensor): internal variables + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: derivative of flow rate with respect to the strain rate + """ + erate, _ = self.flow_rate(s, h, t, T, e) + derate = self.dflow_derate(s, h, t, T, e) + + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] + + return torch.cat( + [ + self.isotropic.dhistory_rate_dtotalrate(s, hiso, t, erate, T, e) + utility.mbmm( + self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e), + derate[..., None, None], + )[...,0], + self.kinematic.dhistory_rate_dtotalrate(s, hkin, t, erate, T, e) + utility.mbmm( + self.kinematic.dhistory_rate_derate(s, hkin, t, erate, T, e), + derate[..., None, None])[...,0], + ], + dim=-1, + ) + + diff --git a/test/test_flowrules.py b/test/test_flowrules.py index 24cc9a6..b39499c 100644 --- a/test/test_flowrules.py +++ b/test/test_flowrules.py @@ -19,7 +19,7 @@ def test_flow_rate(self): nbatch_dim=1, ) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_rate(self): if self.skip: @@ -34,7 +34,7 @@ def test_history_rate(self): nbatch_dim=1, ) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_flow_history(self): if self.skip: @@ -47,7 +47,7 @@ def test_flow_history(self): nbatch_dim=1, ).unsqueeze(1) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_stress(self): if self.skip: @@ -60,7 +60,7 @@ def test_history_stress(self): nbatch_dim=1, ) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_flow_erate(self): exact = self.model.dflow_derate(self.s, self.h, self.t, self.T, self.erate) @@ -70,7 +70,7 @@ def test_flow_erate(self): nbatch_dim=1, ) - self.assertTrue(np.allclose(exact, torch.flatten(numer), rtol=1.0e-2)) + self.assertTrue(np.allclose(exact, torch.flatten(numer), rtol=self.rtol2)) def test_history_erate(self): if self.skip: @@ -83,7 +83,7 @@ def test_history_erate(self): nbatch_dim=1, ) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-2, atol=1e-6)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol2, atol=1e-6)) class CommonFlowRuleBatchBatch: @@ -111,7 +111,7 @@ def test_flow_rate_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_rate_bb(self): if self.skip: @@ -126,7 +126,7 @@ def test_history_rate_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_flow_history_bb(self): if self.skip: @@ -141,7 +141,7 @@ def test_flow_history_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_stress_bb(self): if self.skip: @@ -156,7 +156,7 @@ def test_history_stress_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-4)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_flow_erate_bb(self): @@ -169,7 +169,7 @@ def test_flow_erate_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-2)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol2)) def test_history_erate_bb(self): if self.skip: @@ -184,7 +184,7 @@ def test_history_erate_bb(self): self.assertEqual(exact.shape, numer.shape) - self.assertTrue(np.allclose(exact, numer, rtol=1.0e-2, atol=1e-6)) + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol2, atol=1e-6)) class TestPerfectViscoplasticity( @@ -206,6 +206,29 @@ def setUp(self): self.T = torch.zeros_like(self.t) self.erate = torch.linspace(1e-2, 1e-3, self.nbatch) + self.rtol = 1e-4 + self.rtol2 = 1e-4 + +class TestApproximatePerfectRateIndependentPlasticity( + unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch +): + def setUp(self): + + self.sy = torch.tensor(100.0) + + self.nbatch = 10 + + self.model = flowrules.ApproximatePerfectRateIndependentPlasticity(CP(self.sy)) + + self.s = torch.linspace(90, 100, self.nbatch) + self.h = torch.reshape(torch.linspace(50, 110, self.nbatch), (self.nbatch, 1)) + self.t = torch.ones(self.nbatch) + self.skip = True + self.T = torch.zeros_like(self.t) + self.erate = torch.linspace(1e-2, 1e-3, self.nbatch) + + self.rtol = 1e-4 + self.rtol2 = 1e-4 class TestWrappedRIIsoKinViscoplasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -263,6 +286,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-4 + class TestKocksMeckingRegimeFlowRule( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -332,6 +358,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-2 + class TestSoftKocksMeckingRegimeFlowRule( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -409,6 +438,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-2 + class TestSoftKocksMeckingRegimeFlowRuleComplex( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -485,6 +517,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-2 + class TestIsoKinViscoplasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -527,6 +562,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-4 + def test_kin(self): def fn(i): hp = self.h.clone() @@ -548,6 +586,68 @@ def fn(i): self.assertTrue(np.allclose(i1, i2, rtol=1.0e-4)) +class TestApproximateIsoKinRateIndependentPlasticity( + unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch +): + def setUp(self): + self.sy = torch.tensor(110.0) + + self.nbatch = 10 + + self.R = torch.tensor(101.0) + self.d = torch.tensor(1.3) + self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d)) + + self.C = torch.tensor(1200.0) + self.g = torch.tensor(10.1) + self.kin = hardening.FAKinematicHardeningModel(CP(self.C), CP(self.g)) + + self.model = flowrules.ApproximateIsoKinRateIndependentPlasticity( + CP(self.sy), self.iso, self.kin + ) + + self.s = torch.linspace(380, 140, self.nbatch) + self.h = torch.reshape( + torch.tensor( + np.array( + [ + np.linspace(51, 110, self.nbatch), + np.linspace(-100, 210, self.nbatch)[::-1], + ] + ) + ).T, + (self.nbatch, 2), + ) + + self.t = torch.ones(self.nbatch) + self.T = torch.zeros_like(self.t) + self.erate = torch.linspace(1e-2, 1e-3, self.nbatch) + + self.skip = False + + self.rtol = 1e-3 + self.rtol2 = 1e-2 + + def test_kin(self): + def fn(i): + hp = self.h.clone() + hp[:, 1] = i + return self.model.flow_rate(self.s, hp, self.t, self.T, self.erate)[0] + + i1 = utility.differentiate(fn, self.h[:, 1]) + i2 = self.model.dflow_dkin(self.s, self.h, self.t, self.T, self.erate) + self.assertTrue(np.allclose(i1, i2, rtol=self.rtol)) + + def test_iso(self): + def fn(i): + hp = self.h.clone() + hp[:, 0] = i + return self.model.flow_rate(self.s, hp, self.t, self.T, self.erate)[0] + + i1 = utility.differentiate(fn, self.h[:, 0]) + i2 = self.model.dflow_diso(self.s, self.h, self.t, self.T, self.erate) + + self.assertTrue(np.allclose(i1, i2, rtol=self.rtol)) class TestSuperimposedFlowRate( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -608,6 +708,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-4 + class TestIsoKinChabocheViscoplasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch @@ -651,6 +754,9 @@ def setUp(self): self.skip = False + self.rtol = 1e-4 + self.rtol2 = 1e-4 + def test_kin(self): def fn(i): hp = self.h.clone() From 4cd7f57f481310c6d893fa0a68d68c03b310cbf6 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Mon, 8 Jan 2024 21:19:46 -0600 Subject: [PATCH 08/18] Possible hack fix: --- .../rate_independent_approximations.py | 125 ++++++++ .../stress_strain_control.py | 4 +- pyoptmat/chunktime.py | 9 +- pyoptmat/flowrules.py | 281 +++++++++++------- pyoptmat/models.py | 17 +- pyoptmat/ode.py | 4 + pyoptmat/solvers.py | 4 +- test/test_flowrules.py | 149 +++++----- test/test_models.py | 43 +++ 9 files changed, 446 insertions(+), 190 deletions(-) create mode 100755 examples/structural-material-models/rate_independent_approximations.py diff --git a/examples/structural-material-models/rate_independent_approximations.py b/examples/structural-material-models/rate_independent_approximations.py new file mode 100755 index 0000000..7753dd5 --- /dev/null +++ b/examples/structural-material-models/rate_independent_approximations.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 + +""" + This example demonstrates the ability of pyoptmat to integrate material models under + strain control, stress control, or a combination of both. + + The example first simulates standard, strain-rate controlled tensile curves + for a simple material model for Alloy 617 developed by Messner, Phan, and Sham. These + simulations then take time, strain, and temperature as inputs and outputs the + resulting stresses. + + Then, the example uses the stresses from the strain controlled simulations + to repeat the simulations under stress control. These simulations therefore + take as input time, temperature, and stress and output strain. The plot demonstrates that + the strain controlled and stress controlled simulations produce identical results. + + Finally, the example integrates the model twice for each temperature, once in strain + control and once in stress control, but integrated all the experiments at once. + This example then demos the ability of pyoptmat to integrate both stress and strain + controlled tests in the same vectorized calculation. +""" + +import sys + +sys.path.append("../..") + +import torch + +import matplotlib.pyplot as plt + +from pyoptmat import models, flowrules, temperature, experiments, hardening + +torch.set_default_tensor_type(torch.DoubleTensor) + +if __name__ == "__main__": + # This example illustrates the temperature and rate sensitivity of + # Alloy 617 with the model given in: + # Messner, Phan, and Sham. IJPVP 2019. + + # Setup test conditions + ntemps = 5 + temps = torch.linspace(750, 950, ntemps) + 273.15 + elimits = torch.ones(ntemps) * 0.25 + erates = torch.logspace(-16,-1,ntemps) + + nsteps = 200 + + times, strains, temperatures, cycles = experiments.make_tension_tests( + erates, temps, elimits, nsteps + ) + + # Perfect viscoplastic model + C = torch.tensor(-5.0) + mu = temperature.PolynomialScaling( + torch.tensor([-1.34689305e-02, -5.18806776e00, 7.86708330e04]) + ) + E = temperature.PolynomialScaling( + torch.tensor([-3.48056033e-02, -1.44398964e01, 2.06464967e05]) + ) + sy = temperature.ConstantParameter(200.0) + + C = torch.tensor(12000.0) + g = torch.tensor(10.1) + kh = hardening.FAKinematicHardeningModel(temperature.ConstantParameter(C), temperature.ConstantParameter(g)) + + ih = hardening.Theta0VoceIsotropicHardeningModel(temperature.ConstantParameter(100.0), temperature.ConstantParameter(500.0)) + + fr = flowrules.IsoKinRateIndependentPlasticity( + E, sy, ih, kh + ) + model = models.InelasticModel(E, fr) + integrator = models.ModelIntegrator(model, guess_type = "previous", linesearch = True, block_size = 10) + + stresses = integrator.solve_strain(times, strains, temperatures)[:, :, 0] + + # Now do it backwards! + strains_prime = integrator.solve_stress(times, stresses, temperatures)[:, :, 0] + + plt.figure() + for ei, epi, Ti, si in zip( + strains.T.numpy(), + strains_prime.T.numpy(), + temperatures.T.numpy(), + stresses.T.numpy(), + ): + (l,) = plt.plot(ei, si, label="T = %3.0fK" % Ti[0]) + plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw = 4, alpha = 0.5) + + plt.legend(loc="best") + plt.xlabel("Strain (mm/mm)") + plt.ylabel("Stress (Mpa)") + plt.title("Comparison between strain and stress control") + plt.show() + + # Now do it both ways at once! + big_times = torch.cat((times, times), 1) + big_temperatures = torch.cat((temperatures, temperatures), 1) + big_data = torch.cat((strains, stresses), 1) + big_types = torch.zeros(2 * ntemps) + big_types[ntemps:] = 1 + + both_results = integrator.solve_both( + big_times, big_temperatures, big_data, big_types + ) + + plt.figure() + for i in range(ntemps): + (l,) = plt.plot( + strains[:, i], both_results[:, i, 0], label="T = %3.0fK" % Ti[i] + ) + plt.plot( + both_results[:, i + ntemps, 0], + stresses[:, i], + label=None, + ls="--", + color=l.get_color(), + lw = 4, + alpha = 0.5 + ) + + plt.xlabel("Strain (mm/mm)") + plt.ylabel("Stress (MPa)") + plt.legend(loc="best") + plt.title("Running both strain and stress control at once") + plt.show() diff --git a/examples/structural-material-models/stress_strain_control.py b/examples/structural-material-models/stress_strain_control.py index 051e201..a4234e3 100755 --- a/examples/structural-material-models/stress_strain_control.py +++ b/examples/structural-material-models/stress_strain_control.py @@ -91,7 +91,7 @@ stresses.T.numpy(), ): (l,) = plt.plot(ei, si, label="T = %3.0fK" % Ti[0]) - plt.plot(epi, si, label=None, ls="--", color=l.get_color()) + plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw = 4, alpha = 0.5) plt.legend(loc="best") plt.xlabel("Strain (mm/mm)") @@ -120,7 +120,7 @@ stresses[:, i], label=None, ls="--", - color=l.get_color(), + color=l.get_color(), lw = 4, alpha = 0.5 ) plt.xlabel("Strain (mm/mm)") diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index b7bdb6f..6187bc0 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -58,7 +58,7 @@ def newton_raphson_chunk( R, J = fn(x) nR = torch.norm(R, dim = -1) i += 1 - + if i == miter: if throw_on_fail: raise RuntimeError("Implicit solve did not succeed.") @@ -66,7 +66,7 @@ def newton_raphson_chunk( return x -def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=1e-3, miter = 10): +def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=0.0, miter = 10): """ Backtracking linesearch for the chunk NR algorithm. @@ -93,9 +93,10 @@ def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=1 nR = torch.norm(R, dim = -1)**2.0 crit = nR0**2.0 + 2.0 * c * alpha * torch.einsum('...i,...i', R0, dx) i += 1 - if torch.all(nR <= crit) or i >= miter or torch.all(torch.sqrt(nR) < overall_atol) or torch.all(torch.sqrt(nR)/nR0 < overall_rtol): + if torch.all(nR < crit) or i >= miter or torch.all(torch.sqrt(nR) < overall_atol) or torch.all(torch.sqrt(nR)/nR0 < overall_rtol): break - alpha[nR > crit] /= sigma + alpha[nR >= crit] /= sigma + return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim = -1), alpha diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index a3b2dd3..c9034aa 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -27,7 +27,7 @@ import torch from torch import nn -from pyoptmat import utility +from pyoptmat import utility, solvers class FlowRule(nn.Module): @@ -40,6 +40,9 @@ class FlowRule(nn.Module): def __init__(self): super().__init__() + def setup(self, t, y): + pass + def dflow_dhist(self, s, h, t, T, e): """ The derivative of the flow rate with respect to the internal variables @@ -187,6 +190,10 @@ def __init__( "model must have the same number of internal variables" ) + def setup(self, t, y): + self.model1.setup(t, y) + self.model2.setup(t, y) + @property def nhist(self): """ @@ -452,6 +459,10 @@ def __init__( "model must have the same number of internal variables" ) + def setup(self, t, y): + self.model1.setup(t, y) + self.model2.setup(t, y) + @property def nhist(self): """ @@ -1125,25 +1136,31 @@ def history_rate(self, s, h, t, T, e): """ return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) -class ApproximatePerfectRateIndependentPlasticity(FlowRule): +class PerfectRateIndependentPlasticity(FlowRule): """ - Perfect viscoplasticity defined as - - .. math:: - - \\dot{\\varepsilon}_{in}=\\left(\\frac{\\left|\\sigma\\right|}{\\eta}\\right)^{n}\\operatorname{sign}\\left(\\sigma\\right) - - \\dot{h} = \\emptyset + Perfect rate independent plasticity Args: - n (|TP|): rate sensitivity - eta (|TP|): flow viscosity + sy (|TP|): yield stress """ - def __init__(self, sy, n = 20.0): + def __init__(self, sy): super().__init__() self.sy = sy - self.n = n + + def f(self, s, h, T): + """ + The yield function + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + + Returns: + the value of the yield function... + """ + return torch.abs(s) - self.sy(T) def flow_rate(self, s, h, t, T, e): """ @@ -1162,10 +1179,13 @@ def flow_rate(self, s, h, t, T, e): of the flow rate with respect to stress """ - return ( - torch.abs(e) * utility.macaulay(torch.abs(s) - self.sy(T)) ** self.n * torch.sign(s), - torch.abs(e) * self.n * utility.macaulay(torch.abs(s) - self.sy(T))**(self.n-1) * utility.heaviside(torch.abs(s) - self.sy(T)) * torch.sign(s) - ) + fr = torch.zeros_like(s) + yielding = self.f(s, h, T) > 0 + fr[yielding] = e[yielding] + + dfr = torch.zeros_like(s) + + return fr, dfr def dflow_derate(self, s, h, t, T, e): """ @@ -1185,7 +1205,11 @@ def dflow_derate(self, s, h, t, T, e): torch.tensor: derivative of flow rate with respect to the internal variables """ - return torch.sign(e) * utility.macaulay(torch.abs(s) - self.sy(T)) ** self.n * torch.sign(s) + fr = torch.zeros_like(s) + yielding = self.f(s, h, T) > 0 + fr[yielding] = 1.0 + + return fr @property def nhist(self): @@ -1214,7 +1238,6 @@ def history_rate(self, s, h, t, T, e): """ return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) - class IsoKinViscoplasticity(FlowRule): """ Viscoplasticity with isotropic and kinematic hardening, defined as @@ -1516,7 +1539,7 @@ def dhist_derate(self, s, h, t, T, e): ) -class ApproximateIsoKinRateIndependentPlasticity(FlowRule): +class IsoKinRateIndependentPlasticity(FlowRule): """ Viscoplasticity with isotropic and kinematic hardening, defined as @@ -1546,70 +1569,111 @@ class ApproximateIsoKinRateIndependentPlasticity(FlowRule): kinematic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the kinematic hardening model """ - def __init__(self, sy, isotropic, kinematic, n = 20.0): + def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10): super().__init__() + self. E = E self.isotropic = isotropic self.kinematic = kinematic - self.n = n self.sy = sy + self.soffset = soffset - def flow_rate(self, s, h, t, T, e): + self.yielding = None + + def setup(self, t, y): + self.yielding = None + + def update_yielding(self, s, h, T): + if self.yielding is None or self.yielding.shape != s.shape: + self.yielding = self.f(s, h, T) > 0 + else: + self.yielding = torch.logical_or(self.yielding,self.f(s, h, T) > 0) + + def f(self, s, h, T): """ - The flow rate itself and the derivative with respect to stress + The yield function Args: s (torch.tensor): stress h (torch.tensor): history - t (torch.tensor): time T (torch.tensor): temperature - e (torch.tensor): total strain rate Returns: - tuple(torch.tensor, torch.tensor): the flow rate and the derivative of the flow rate with - respect to stress + the value of the yield function... """ ih = self.isotropic.value(h[..., : self.isotropic.nhist]) kh = self.kinematic.value(h[..., self.isotropic.nhist :]) + return torch.abs(s-kh) - self.sy(T) - ih - return ( - torch.abs(e) * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) - ** self.n - * torch.sign(s - kh), - torch.abs(e) * self.n - * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) - ** (self.n - 1), - ) + def ep_residual(self, ep, s, h, t, T, e): + """ + The residual function to solve for the consistency parameter, here just + expanded to the plastic flow rate as it's a scalar. + """ + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] - def dflow_diso(self, s, h, t, T, e): + kh = self.kinematic.value(hkin) + + i_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.history_rate(s, hiso, t, ep, T, e).unsqueeze(-1)).squeeze(-1).squeeze(-1) + k_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.history_rate(s, hkin, t, ep, T, e).unsqueeze(-1)).squeeze(-1).squeeze(-1) + + return torch.sign(s - kh + self.soffset) * self.E(T) * (e - ep) - torch.sign(s - kh + self.soffset) * k_dot - i_dot + + def ep_jacobian_ep(self, ep, s, h, t, T, e): """ - The derivative of the flow rate with respect to the isotropic hardening + The Jacobian for the plastic flow rate residual with respect to the plastic strain rate + """ + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] - Args: - s (torch.tensor): stress - h (torch.tensor): history - t (torch.tensor): time - T (torch.tensor): temperature - e (torch.tensor): total strain rate + kh = self.kinematic.value(hkin) + + di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_derate(s, hiso, t, ep, T, e)).squeeze(-1) + dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_derate(s, hkin, t, ep, T, e)).squeeze(-1) - Returns: - torch.tensor: the derivative of the flow rate + return -(torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) - torch.sign(s - kh + self.soffset).unsqueeze(-1) * dk_dot - di_dot + + def ep_jacobian_s(self, ep, s, h, t, T, e): """ - ih = self.isotropic.value(h[..., : self.isotropic.nhist]) - kh = self.kinematic.value( - h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] - ) + The Jacobian for the plastic flow rate residual with respect to the stress + """ + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] - iv = (torch.abs(s - kh) - self.sy(T) - ih) + kh = self.kinematic.value(hkin) + + di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_dstress(s, hiso, t, ep, T, e).unsqueeze(-1)).squeeze(-1) + dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_dstress(s, hkin, t, ep, T, e).unsqueeze(-1)).squeeze(-1) - return ( - -torch.abs(e) * self.n - * utility.macaulay(iv) ** (self.n - 1) - * torch.sign(s - kh) - ) + return -torch.sign(s - kh + self.soffset).unsqueeze(-1) * dk_dot - di_dot - def dflow_dkin(self, s, h, t, T, e): + def ep_jacobian_h(self, ep, s, h, t, T, e): """ - The derivative of the flow rate with respect to the kinematic hardening + The Jacobian for the plastic flow rate residual with respect to the stress + """ + hiso = h[..., : self.isotropic.nhist] + hkin = h[..., self.isotropic.nhist :] + + kh = self.kinematic.value(hkin) + + di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_dhistory(s, hiso, t, ep, T, e)).squeeze(-2) + dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_dhistory(s, hkin, t, ep, T, e)).squeeze(-2) + + return torch.cat([-di_dot, -torch.sign(s - kh).unsqueeze(-1) * dk_dot], dim = -1) + + def ep_jacobian_e(self, ep, s, h, t, T, e): + """ + The Jacobian for the plastic flow rate residual with respect to the total strain rate + """ + hkin = h[..., self.isotropic.nhist :] + + kh = self.kinematic.value(hkin) + + return (torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) + + def flow_rate(self, s, h, t, T, e): + """ + The flow rate itself and the derivative with respect to stress Args: s (torch.tensor): stress @@ -1619,18 +1683,22 @@ def dflow_dkin(self, s, h, t, T, e): e (torch.tensor): total strain rate Returns: - torch.tensor: the derivative of the flow rate + tuple(torch.tensor, torch.tensor): the flow rate and the derivative of the flow rate with + respect to stress """ - ih = self.isotropic.value(h[..., : self.isotropic.nhist]) - kh = self.kinematic.value( - h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] - ) + fr = torch.zeros_like(e) + self.update_yielding(s, h, T) + + def RJ(x): + return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1) - return ( - -torch.abs(e) * self.n - * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih)) - ** (self.n - 1) - ) + pfr = solvers.scalar_newton(RJ, torch.zeros_like(e)) + fr[self.yielding] = pfr[self.yielding] + dfr = torch.zeros_like(fr) + dfr[self.yielding] = -(self.ep_jacobian_s(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding] + + return fr, dfr + def dflow_derate(self, s, h, t, T, e): """ @@ -1650,12 +1718,32 @@ def dflow_derate(self, s, h, t, T, e): torch.tensor: derivative of flow rate with respect to the internal variables """ - ih = self.isotropic.value(h[..., : self.isotropic.nhist]) - kh = self.kinematic.value( - h[..., self.isotropic.nhist : self.isotropic.nhist + self.kinematic.nhist] - ) + pfr = self.flow_rate(s, h, t, T, e)[0] + dfr = torch.zeros_like(e) + dfr[self.yielding] = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding] + + return dfr + + def dflow_dhist(self, s, h, t, T, e): + """ + The derivative of the flow rate with respect to the internal variables + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + + Returns: + torch.tensor: the derivative of the flow rate + """ + pfr = self.flow_rate(s, h, t, T, e)[0] + dfr = torch.zeros_like(h) + + dfr[self.yielding] = -(self.ep_jacobian_h(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e))[self.yielding] - return torch.sign(e) * utility.macaulay((torch.abs(s - kh) - self.sy(T) - ih))**self.n * torch.sign(s - kh) + return dfr.unsqueeze(-2) @property def nhist(self): @@ -1699,6 +1787,10 @@ def history_rate(self, s, h, t, T, e): dim=-1, ) + df_dh = self.dflow_dhist(s, h, t, T, e) + df_di = df_dh[...,:self.isotropic.nhist] + df_dk = df_dh[...,self.isotropic.nhist:] + # Jacobian contribution hdiv = torch.cat( [ @@ -1709,14 +1801,10 @@ def history_rate(self, s, h, t, T, e): self.isotropic.dhistory_rate_derate( s, hiso, t, erate, T, e ), - utility.mbmm( - self.dflow_diso(s, h, t, T, e)[..., None, None], - self.isotropic.dvalue(hiso)[..., None, :], - ), + df_di ), self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e) - * self.dflow_dkin(s, h, t, T, e)[..., None, None] - * self.kinematic.dvalue(hkin)[..., None, :], + * df_dk, ], dim=-1, ), @@ -1726,20 +1814,14 @@ def history_rate(self, s, h, t, T, e): self.kinematic.dhistory_rate_derate( s, hkin, t, erate, T, e ), - utility.mbmm( - self.dflow_diso(s, h, t, T, e)[..., None, None], - self.isotropic.dvalue(hiso)[..., None, :], - ), + df_di ), self.kinematic.dhistory_rate_dhistory(s, hkin, t, erate, T, e) + utility.mbmm( self.kinematic.dhistory_rate_derate( s, hkin, t, erate, T, e ), - utility.mbmm( - self.dflow_dkin(s, h, t, T, e)[..., None, None], - self.kinematic.dvalue(hkin)[..., None, :], - ), + df_dk, ), ], dim=-1, @@ -1750,31 +1832,6 @@ def history_rate(self, s, h, t, T, e): return hrate, hdiv - def dflow_dhist(self, s, h, t, T, e): - """ - The derivative of the flow rate with respect to the internal variables - - Args: - s (torch.tensor): stress - h (torch.tensor): history - t (torch.tensor): time - T (torch.tensor): temperature - e (torch.tensor): total strain rate - - Returns: - torch.tensor: the derivative of the flow rate - """ - hiso = h[..., : self.isotropic.nhist] - hkin = h[..., self.isotropic.nhist :] - - return torch.cat( - [ - self.dflow_diso(s, h, t, T, e)[..., None] * self.isotropic.dvalue(hiso), - self.dflow_dkin(s, h, t, T, e)[..., None] * self.kinematic.dvalue(hkin), - ], - dim=-1, - ).unsqueeze(-2) - def dhist_dstress(self, s, h, t, T, e): """ The derivative of the history rate with respect to the stress diff --git a/pyoptmat/models.py b/pyoptmat/models.py index a56054e..93cdfe1 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -72,6 +72,9 @@ def nhist(self): """ return 1 + self.flowrule.nhist + def setup(self, t, y): + self.flowrule.setup(t, y) + def forward(self, t, y, erate, T): """ Return the rate equations for the strain-based version of the model @@ -133,7 +136,7 @@ def forward(self, t, y, erate, T): # Logically we should return the derivative wrt T, but right now # we're not going to use it Trate = torch.zeros_like(y) - + return result, dresult, drate, Trate @@ -452,6 +455,12 @@ def __init__(self, model, times, rates, base, temps, control, bisect_first = Fal bisect_first = bisect_first ) + def setup(self, t, y): + if torch.any(self.econtrol): + self.emodel.setup(t[...,self.econtrol], y[...,self.econtrol,:]) + if torch.any(self.scontrol): + self.smodel.setup(t[...,self.scontrol], y[...,self.scontrol,:]) + def forward(self, t, y): """ Evaluate both strain and stress control and paste into the right @@ -496,6 +505,9 @@ def __init__(self, model, erate_fn, T_fn, *args, **kwargs): self.erate_fn = erate_fn self.T_fn = T_fn + def setup(self, t, y): + self.model.setup(t, y) + def forward(self, t, y): """ Strain rate as a function of t and state @@ -529,6 +541,9 @@ def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate self.max_erate = max_erate self.bisect_first = bisect_first + def setup(self, t, y): + self.model.setup(t, y) + def forward(self, t, y): """ Stress rate as a function of t and state diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index 007de21..716cf27 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -520,6 +520,10 @@ def block_update(self, t, t_start, y_start, func, y_guess): """ # Various useful sizes n = t.shape[0] # Number of time steps to do at once + + self.func.setup(torch.vstack((t_start.unsqueeze(0), t)), + torch.vstack((torch.zeros_like(y_start).unsqueeze(0), y_guess.reshape(self.batch_size, n, self.prob_size).transpose(0, 1))) + + y_start.unsqueeze(0).expand(n + 1, -1, -1)) def RJ(dy): # Make things into a more rational shape diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index a6014c2..bcaf7d4 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -61,9 +61,9 @@ def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100): for i in range(miter): if torch.all(torch.abs(R) < atol): break - + x -= R / J - + R, J = fn(x) else: warnings.warn("Scalar implicit solve did not succeed. Results may be inaccurate...") diff --git a/test/test_flowrules.py b/test/test_flowrules.py index b39499c..e9a9d9b 100644 --- a/test/test_flowrules.py +++ b/test/test_flowrules.py @@ -18,7 +18,7 @@ def test_flow_rate(self): self.s, nbatch_dim=1, ) - + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_rate(self): @@ -46,7 +46,7 @@ def test_flow_history(self): self.h, nbatch_dim=1, ).unsqueeze(1) - + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_stress(self): @@ -149,7 +149,7 @@ def test_history_stress_bb(self): s, h, t, T, erate = self.expand_all() - exact = self.model.dhist_dstress(s, h, t, T, self.erate) + exact = self.model.dhist_dstress(s, h, t, T, erate) numer = utility.batch_differentiate( lambda x: self.model.history_rate(x, h, t, T, erate)[0], s, nbatch_dim=2 ) @@ -209,7 +209,7 @@ def setUp(self): self.rtol = 1e-4 self.rtol2 = 1e-4 -class TestApproximatePerfectRateIndependentPlasticity( +class TestPerfectRateIndependentPlasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch ): def setUp(self): @@ -218,9 +218,9 @@ def setUp(self): self.nbatch = 10 - self.model = flowrules.ApproximatePerfectRateIndependentPlasticity(CP(self.sy)) + self.model = flowrules.PerfectRateIndependentPlasticity(CP(self.sy)) - self.s = torch.linspace(90, 100, self.nbatch) + self.s = torch.linspace(91, 121, self.nbatch) self.h = torch.reshape(torch.linspace(50, 110, self.nbatch), (self.nbatch, 1)) self.t = torch.ones(self.nbatch) self.skip = True @@ -586,69 +586,6 @@ def fn(i): self.assertTrue(np.allclose(i1, i2, rtol=1.0e-4)) -class TestApproximateIsoKinRateIndependentPlasticity( - unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch -): - def setUp(self): - self.sy = torch.tensor(110.0) - - self.nbatch = 10 - - self.R = torch.tensor(101.0) - self.d = torch.tensor(1.3) - self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d)) - - self.C = torch.tensor(1200.0) - self.g = torch.tensor(10.1) - self.kin = hardening.FAKinematicHardeningModel(CP(self.C), CP(self.g)) - - self.model = flowrules.ApproximateIsoKinRateIndependentPlasticity( - CP(self.sy), self.iso, self.kin - ) - - self.s = torch.linspace(380, 140, self.nbatch) - self.h = torch.reshape( - torch.tensor( - np.array( - [ - np.linspace(51, 110, self.nbatch), - np.linspace(-100, 210, self.nbatch)[::-1], - ] - ) - ).T, - (self.nbatch, 2), - ) - - self.t = torch.ones(self.nbatch) - self.T = torch.zeros_like(self.t) - self.erate = torch.linspace(1e-2, 1e-3, self.nbatch) - - self.skip = False - - self.rtol = 1e-3 - self.rtol2 = 1e-2 - - def test_kin(self): - def fn(i): - hp = self.h.clone() - hp[:, 1] = i - return self.model.flow_rate(self.s, hp, self.t, self.T, self.erate)[0] - - i1 = utility.differentiate(fn, self.h[:, 1]) - i2 = self.model.dflow_dkin(self.s, self.h, self.t, self.T, self.erate) - self.assertTrue(np.allclose(i1, i2, rtol=self.rtol)) - - def test_iso(self): - def fn(i): - hp = self.h.clone() - hp[:, 0] = i - return self.model.flow_rate(self.s, hp, self.t, self.T, self.erate)[0] - - i1 = utility.differentiate(fn, self.h[:, 0]) - i2 = self.model.dflow_diso(self.s, self.h, self.t, self.T, self.erate) - - self.assertTrue(np.allclose(i1, i2, rtol=self.rtol)) - class TestSuperimposedFlowRate( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch ): @@ -777,3 +714,77 @@ def fn(i): i2 = self.model.dflow_diso(self.s, self.h, self.t, self.T, self.erate) self.assertTrue(np.allclose(i1, i2, rtol=1.0e-4)) + + +class TestIsoKinChabocheRIPlasticity( + unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch +): + def setUp(self): + self.sy = torch.tensor(110.0) + self.E = torch.tensor(50000.0) + + self.nbatch = 10 + + self.R = torch.tensor(101.0) + self.d = torch.tensor(1.3) + self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d)) + + self.C = torch.tensor([100.0, 1000, 1500]) + self.g = torch.tensor([1.2, 100, 50]) + self.kin = hardening.ChabocheHardeningModel(CP(self.C), CP(self.g)) + + self.model = flowrules.IsoKinRateIndependentPlasticity( + CP(self.E), CP(self.sy), self.iso, self.kin + ) + + self.s = torch.linspace(160, 240, self.nbatch) + self.h = torch.reshape( + torch.tensor( + np.array( + [ + np.linspace(51, 110, self.nbatch), + np.linspace(-10, 21, self.nbatch)[::-1], + np.linspace(0, 2, self.nbatch), + np.linspace(-2, 0, self.nbatch), + ] + ) + ).T, + (self.nbatch, 4), + ) + self.t = torch.ones(self.nbatch) + self.T = torch.zeros_like(self.t) + self.erate = torch.linspace(1e-2, 1e-3, self.nbatch) + + self.skip = False + + self.rtol = 1e-4 + self.rtol2 = 1e-4 + + self.ep_guess = self.erate.clone() + + def test_d_residual_dep(self): + exact = self.model.ep_jacobian_ep(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) + numer = utility.batch_differentiate(lambda x: self.model.ep_residual(x, self.s, self.h, self.t, self.T, self.erate), self.ep_guess).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_d_residual_ds(self): + exact = self.model.ep_jacobian_s(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) + numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, x, self.h, self.t, self.T, self.erate), self.s).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_d_residual_dh(self): + exact = self.model.ep_jacobian_h(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) + numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, x, self.t, self.T, self.erate), self.h) + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_d_residual_de(self): + exact = self.model.ep_jacobian_e(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) + numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, self.h, self.t, self.T, x), self.erate).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_definiition_flow_rate(self): + vals = self.model.flow_rate(self.s, self.h, self.t, self.T, self.erate)[0] + yielding = self.model.f(self.s, self.h, self.T) > 0 + r = self.model.ep_residual(vals, self.s, self.h, self.t, self.T, self.erate) + self.assertTrue(torch.allclose(r[yielding], torch.tensor(0.0))) + diff --git a/test/test_models.py b/test/test_models.py index 5cd51ef..61e87a0 100644 --- a/test/test_models.py +++ b/test/test_models.py @@ -327,6 +327,48 @@ def setUp(self): self.t = self.times[2] self.step = 2 +class TestIsoKinRIPlasticity(unittest.TestCase, CommonModel, CommonModelBatchBatch): + def setUp(self): + self.E = torch.tensor(100000.0) + self.sy = torch.tensor(10.0) + + self.R = torch.tensor(101.0) + self.d = torch.tensor(1.3) + self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d)) + + self.C = torch.tensor(12000.0) + self.g = torch.tensor(10.1) + self.kin = hardening.FAKinematicHardeningModel(CP(self.C), CP(self.g)) + + self.flowrule = flowrules.IsoKinRateIndependentPlasticity( + CP(self.E), CP(self.sy), self.iso, self.kin + ) + self.model = models.InelasticModel(CP(self.E), self.flowrule) + + self.times = torch.transpose( + torch.tensor(np.array([np.linspace(0, 1, 4) for i in range(3)])), 1, 0 + ) + self.strains = torch.transpose( + torch.tensor(np.array([np.linspace(0, 1, 4) for i in range(3)])), 1, 0 + ) + self.temperatures = torch.zeros_like(self.times) + self.stresses = ( + torch.transpose( + torch.tensor(np.array([np.linspace(0, 1, 4) for i in range(3)])), 1, 0 + ) + * 200 + ) + + self.state_strain = ( + torch.tensor([[90.0, 30.0, 10.0], [100.0, 10.0, 15.0], [101.0, 50.0, 60.0]]) + ) / 3.0 + self.state_stress = ( + torch.tensor([[0.05, 30.0, 10.0], [0.07, 10.0, 15.0], [0.08, 50.0, 60.0]]) + ) / 3.0 + + self.t = self.times[2] + self.step = 2 + class TestIsoKinViscoplasticityRecovery( unittest.TestCase, CommonModel, CommonModelBatchBatch @@ -492,3 +534,4 @@ def setUp(self): self.t = self.times[2] self.step = 2 + From c8f87d709cf9740927a133fbb954dd7ddb8474f3 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 11 Jan 2024 18:08:29 -0600 Subject: [PATCH 09/18] Sooo close --- pyoptmat/chunktime.py | 1 + pyoptmat/flowrules.py | 251 ++++++++++++++++++++++++++++++++--------- pyoptmat/models.py | 21 +--- pyoptmat/ode.py | 4 - pyoptmat/solvers.py | 4 +- test/test_flowrules.py | 40 +++++-- test/test_models.py | 16 ++- 7 files changed, 250 insertions(+), 87 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index 6187bc0..f9f3ec4 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -50,6 +50,7 @@ def newton_raphson_chunk( i = 0 while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))): + print(i, torch.max(nR)) dx = solver.solve(J, R) if linesearch: x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R, rtol, atol) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index c9034aa..595d615 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -40,9 +40,6 @@ class FlowRule(nn.Module): def __init__(self): super().__init__() - def setup(self, t, y): - pass - def dflow_dhist(self, s, h, t, T, e): """ The derivative of the flow rate with respect to the internal variables @@ -190,10 +187,6 @@ def __init__( "model must have the same number of internal variables" ) - def setup(self, t, y): - self.model1.setup(t, y) - self.model2.setup(t, y) - @property def nhist(self): """ @@ -459,10 +452,6 @@ def __init__( "model must have the same number of internal variables" ) - def setup(self, t, y): - self.model1.setup(t, y) - self.model2.setup(t, y) - @property def nhist(self): """ @@ -1541,15 +1530,25 @@ def dhist_derate(self, s, h, t, T, e): class IsoKinRateIndependentPlasticity(FlowRule): """ - Viscoplasticity with isotropic and kinematic hardening, defined as + An approximation to rate independent plasticity. The model is defined by + + .. math:: + + \\dot{\\varepsilon}_{in} = \\xi(f) \\dot{\\varepsilon}_{p,ri} + + where :math:`\\xi` is a sigmoid function, :math:`f` is a flow surface + and :math:`\\dot{\\varepsilon}_{p,ri}` is the rate independent plastic flow + rate, as defined by the classical consistency conditions. + + This function uses the flow surface .. math:: - \\dot{\\varepsilon}_{in}=\\left\\langle \\frac{\\left|\\sigma-x\\right|-s_{0}-k}{\\eta}\\right\\rangle ^{n}\\operatorname{sign}\\left(\\sigma-X\\right) + f = \\left|\\sigma - x\\right| - \\sigma_y - k and where the :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects determine the - history rate. + history rate of isotropic (:math:`k`) and kinematic (:math:`x`) hardening. The :py:class:`pyoptmat.hardening.IsotropicHardeningModel` and :py:class:`pyoptmat.hardening.KinematicHardeningModel` objects each define both a @@ -1562,14 +1561,17 @@ class IsoKinRateIndependentPlasticity(FlowRule): the information this class needs to provide. Args: - n (|TP|): rate sensitivity - eta (|TP|): flow viscosity - s0 (|TP|): initial value of flow stress (i.e. the threshold stress) + E (|TP|): young's modulus, needed to define the rate independent flow rate + sy (|TP|): yield stress isotropic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the isotropic hardening model kinematic (:py:class:`pyoptmat.hardening.IsotropicHardeningModel`): object providing the kinematic hardening model + + Keyword Args: + soffset (float): small offset to the stress in the yield surface to avoid a singularity at zero + s (float): scale factor for the sigmoid function, controls the amount of smoothing at the onset of plasticity """ - def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10): + def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10, s = 1.0): super().__init__() self. E = E self.isotropic = isotropic @@ -1577,16 +1579,7 @@ def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10): self.sy = sy self.soffset = soffset - self.yielding = None - - def setup(self, t, y): - self.yielding = None - - def update_yielding(self, s, h, T): - if self.yielding is None or self.yielding.shape != s.shape: - self.yielding = self.f(s, h, T) > 0 - else: - self.yielding = torch.logical_or(self.yielding,self.f(s, h, T) > 0) + self.s = s def f(self, s, h, T): """ @@ -1598,16 +1591,60 @@ def f(self, s, h, T): T (torch.tensor): temperature Returns: - the value of the yield function... + the value of the yield function """ ih = self.isotropic.value(h[..., : self.isotropic.nhist]) kh = self.kinematic.value(h[..., self.isotropic.nhist :]) return torch.abs(s-kh) - self.sy(T) - ih + def df_ds(self, s, h, T): + """ + The derivative of the yield function with respect to stress + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + + Returns: + the derivative value + """ + kh = self.kinematic.value(h[..., self.isotropic.nhist :]) + return torch.sign(s-kh) + + def df_dh(self, s, h, T): + """ + The derivative of the yield function with respect to the history + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + + Returns: + the derivative value + """ + kh = self.kinematic.value(h[..., self.isotropic.nhist :]) + di = -self.isotropic.dvalue(h[..., : self.isotropic.nhist]) + dk = -torch.sign(s-kh).unsqueeze(-1) * self.kinematic.dvalue(h[..., self.isotropic.nhist :]) + + return torch.cat([di,dk], axis = -1) + def ep_residual(self, ep, s, h, t, T, e): """ The residual function to solve for the consistency parameter, here just expanded to the plastic flow rate as it's a scalar. + + Args: + ep (torch.tensor): current values of the plastic flow rate + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): current time + T (torch.tensor): temperature + e (torch.tensor): current total strain rate + + Returns: + the residual value """ hiso = h[..., : self.isotropic.nhist] hkin = h[..., self.isotropic.nhist :] @@ -1621,7 +1658,18 @@ def ep_residual(self, ep, s, h, t, T, e): def ep_jacobian_ep(self, ep, s, h, t, T, e): """ - The Jacobian for the plastic flow rate residual with respect to the plastic strain rate + The Jacobian of the consistency residual with respect to the plastic strain rate. + + Args: + ep (torch.tensor): current values of the plastic flow rate + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): current time + T (torch.tensor): temperature + e (torch.tensor): current total strain rate + + Returns: + the derivative value """ hiso = h[..., : self.isotropic.nhist] hkin = h[..., self.isotropic.nhist :] @@ -1635,7 +1683,18 @@ def ep_jacobian_ep(self, ep, s, h, t, T, e): def ep_jacobian_s(self, ep, s, h, t, T, e): """ - The Jacobian for the plastic flow rate residual with respect to the stress + The Jacobian of the consistency residual with respect to the stress + + Args: + ep (torch.tensor): current values of the plastic flow rate + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): current time + T (torch.tensor): temperature + e (torch.tensor): current total strain rate + + Returns: + the derivative value """ hiso = h[..., : self.isotropic.nhist] hkin = h[..., self.isotropic.nhist :] @@ -1649,7 +1708,18 @@ def ep_jacobian_s(self, ep, s, h, t, T, e): def ep_jacobian_h(self, ep, s, h, t, T, e): """ - The Jacobian for the plastic flow rate residual with respect to the stress + The Jacobian of the consistency residual with respect to the history + + Args: + ep (torch.tensor): current values of the plastic flow rate + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): current time + T (torch.tensor): temperature + e (torch.tensor): current total strain rate + + Returns: + the derivative value """ hiso = h[..., : self.isotropic.nhist] hkin = h[..., self.isotropic.nhist :] @@ -1663,7 +1733,18 @@ def ep_jacobian_h(self, ep, s, h, t, T, e): def ep_jacobian_e(self, ep, s, h, t, T, e): """ - The Jacobian for the plastic flow rate residual with respect to the total strain rate + The Jacobian of the consistency residual with respect to the total strain rate + + Args: + ep (torch.tensor): current values of the plastic flow rate + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): current time + T (torch.tensor): temperature + e (torch.tensor): current total strain rate + + Returns: + the derivative value """ hkin = h[..., self.isotropic.nhist :] @@ -1671,6 +1752,79 @@ def ep_jacobian_e(self, ep, s, h, t, T, e): return (torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) + def sig(self, x): + """ + Our chosen sigmoid function + + .. math:: + + \\xi = \\frac{\\tanh(s x) + 1}{2} + + with :math:`s` the scale factor + + Args: + x (torch.tensor): input to sigmoid + """ + return (torch.tanh(self.s * x) + 1.0)/2.0 + + def dsig(self, x): + """ + Derivative of the sigmoid + + Args: + x (torch.tensor): input to sigmoid + """ + return 0.5*self.s / torch.cosh(self.s * x)**2.0 + + def mix_fn(self, s, h, T): + """ + The mixture function :math:`\\xi(f)` + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + """ + return self.sig(self.f(s, h, T)) + + def dmix_fn_ds(self, s, h, T): + """ + The derivative of the mixture function with respect to the stress + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + """ + return self.dsig(self.f(s, h, T)) * self.df_ds(s, h, T) + + def dmix_fn_dh(self, s, h, T): + """ + The derivative of the mixture function with respect to the history + + Args: + s (torch.tensor): stress + h (torch.tensor): history + T (torch.tensor): temperature + """ + return self.dsig(self.f(s, h, T)).unsqueeze(-1) * self.df_dh(s, h, T) + + def plastic_rate(self, s, h, t, T, e): + """ + Solve for the plastic strain rate that meets the consistency criteria + + Args: + s (torch.tensor): stress + h (torch.tensor): history + t (torch.tensor): time + T (torch.tensor): temperature + e (torch.tensor): total strain rate + """ + def RJ(x): + return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1) + + return solvers.scalar_newton(RJ, torch.zeros_like(e)) + def flow_rate(self, s, h, t, T, e): """ The flow rate itself and the derivative with respect to stress @@ -1686,17 +1840,13 @@ def flow_rate(self, s, h, t, T, e): tuple(torch.tensor, torch.tensor): the flow rate and the derivative of the flow rate with respect to stress """ - fr = torch.zeros_like(e) - self.update_yielding(s, h, T) - - def RJ(x): - return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1) - - pfr = solvers.scalar_newton(RJ, torch.zeros_like(e)) - fr[self.yielding] = pfr[self.yielding] - dfr = torch.zeros_like(fr) - dfr[self.yielding] = -(self.ep_jacobian_s(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding] + # Solve for the plastic flow rate + mf = self.mix_fn(s, h, T) + bfr = self.plastic_rate(s, h, t, T, e) + fr = mf * bfr + dfr = -(self.ep_jacobian_s(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)).squeeze(-1) * mf + self.dmix_fn_ds(s, h, T) * bfr + return fr, dfr @@ -1718,9 +1868,9 @@ def dflow_derate(self, s, h, t, T, e): torch.tensor: derivative of flow rate with respect to the internal variables """ - pfr = self.flow_rate(s, h, t, T, e)[0] - dfr = torch.zeros_like(e) - dfr[self.yielding] = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1)[self.yielding] + pfr = self.plastic_rate(s, h, t, T, e) + mf = self.mix_fn(s, h, T) + dfr = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1) * mf return dfr @@ -1738,11 +1888,10 @@ def dflow_dhist(self, s, h, t, T, e): Returns: torch.tensor: the derivative of the flow rate """ - pfr = self.flow_rate(s, h, t, T, e)[0] - dfr = torch.zeros_like(h) - - dfr[self.yielding] = -(self.ep_jacobian_h(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e))[self.yielding] - + bfr = self.plastic_rate(s, h, t, T, e) + mf = self.mix_fn(s, h, T) + dfr = -(self.ep_jacobian_h(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)) * mf.unsqueeze(-1) + self.dmix_fn_dh(s, h, T) * bfr.unsqueeze(-1) + return dfr.unsqueeze(-2) @property diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 93cdfe1..09f6660 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -72,9 +72,6 @@ def nhist(self): """ return 1 + self.flowrule.nhist - def setup(self, t, y): - self.flowrule.setup(t, y) - def forward(self, t, y, erate, T): """ Return the rate equations for the strain-based version of the model @@ -91,8 +88,8 @@ def forward(self, t, y, erate, T): d_y_dot_d_erate:(nbatch,1+nhist) Jacobian wrt the strain rate d_y_dot_d_T: (nbatch,1+nhist) derivative wrt temperature (unused) """ - stress = y[..., 0].clone() - h = y[..., 1:].clone() + stress = y[..., 0] + h = y[..., 1:] frate, dfrate = self.flowrule.flow_rate(stress, h, t, T, erate) hrate, dhrate = self.flowrule.history_rate(stress, h, t, T, erate) @@ -455,12 +452,6 @@ def __init__(self, model, times, rates, base, temps, control, bisect_first = Fal bisect_first = bisect_first ) - def setup(self, t, y): - if torch.any(self.econtrol): - self.emodel.setup(t[...,self.econtrol], y[...,self.econtrol,:]) - if torch.any(self.scontrol): - self.smodel.setup(t[...,self.scontrol], y[...,self.scontrol,:]) - def forward(self, t, y): """ Evaluate both strain and stress control and paste into the right @@ -505,9 +496,6 @@ def __init__(self, model, erate_fn, T_fn, *args, **kwargs): self.erate_fn = erate_fn self.T_fn = T_fn - def setup(self, t, y): - self.model.setup(t, y) - def forward(self, t, y): """ Strain rate as a function of t and state @@ -541,9 +529,6 @@ def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate self.max_erate = max_erate self.bisect_first = bisect_first - def setup(self, t, y): - self.model.setup(t, y) - def forward(self, t, y): """ Stress rate as a function of t and state @@ -559,7 +544,7 @@ def forward(self, t, y): def RJ(erate): ydot, _, Je, _ = self.model(t, torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), - erate.detach(), cT) + erate, cT) R = ydot[..., 0] - csr J = Je[..., 0] diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index 716cf27..ba8940b 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -521,10 +521,6 @@ def block_update(self, t, t_start, y_start, func, y_guess): # Various useful sizes n = t.shape[0] # Number of time steps to do at once - self.func.setup(torch.vstack((t_start.unsqueeze(0), t)), - torch.vstack((torch.zeros_like(y_start).unsqueeze(0), y_guess.reshape(self.batch_size, n, self.prob_size).transpose(0, 1))) - + y_start.unsqueeze(0).expand(n + 1, -1, -1)) - def RJ(dy): # Make things into a more rational shape dy = dy.reshape(self.batch_size, n, self.prob_size).transpose(0, 1) diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index bcaf7d4..8deed4a 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -62,12 +62,12 @@ def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100): if torch.all(torch.abs(R) < atol): break - x -= R / J + x = x - R / J R, J = fn(x) else: warnings.warn("Scalar implicit solve did not succeed. Results may be inaccurate...") - + return x def scalar_bisection_newton(fn, a, b, atol = 1.0e-6, miter = 100, biter = 10): diff --git a/test/test_flowrules.py b/test/test_flowrules.py index e9a9d9b..65c54dd 100644 --- a/test/test_flowrules.py +++ b/test/test_flowrules.py @@ -18,7 +18,7 @@ def test_flow_rate(self): self.s, nbatch_dim=1, ) - + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_rate(self): @@ -46,7 +46,7 @@ def test_flow_history(self): self.h, nbatch_dim=1, ).unsqueeze(1) - + self.assertTrue(np.allclose(exact, numer, rtol=self.rtol)) def test_history_stress(self): @@ -734,7 +734,7 @@ def setUp(self): self.kin = hardening.ChabocheHardeningModel(CP(self.C), CP(self.g)) self.model = flowrules.IsoKinRateIndependentPlasticity( - CP(self.E), CP(self.sy), self.iso, self.kin + CP(self.E), CP(self.sy), self.iso, self.kin, s = 0.01 ) self.s = torch.linspace(160, 240, self.nbatch) @@ -782,9 +782,33 @@ def test_d_residual_de(self): numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, self.h, self.t, self.T, x), self.erate).unsqueeze(-1) self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) - def test_definiition_flow_rate(self): - vals = self.model.flow_rate(self.s, self.h, self.t, self.T, self.erate)[0] - yielding = self.model.f(self.s, self.h, self.T) > 0 - r = self.model.ep_residual(vals, self.s, self.h, self.t, self.T, self.erate) - self.assertTrue(torch.allclose(r[yielding], torch.tensor(0.0))) + def test_d_sig(self): + x = torch.linspace(-1, 1, 10) + exact = self.model.dsig(x) + numer = utility.batch_differentiate(lambda x: self.model.sig(x), x) + + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_df_ds(self): + exact = self.model.df_ds(self.s, self.h, self.T) + numer = utility.batch_differentiate(lambda x: self.model.f(x, self.h, self.T), self.s) + + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_dmix_fn_ds(self): + exact = self.model.dmix_fn_ds(self.s, self.h, self.T) + numer = utility.batch_differentiate(lambda x: self.model.mix_fn(x, self.h, self.T), self.s) + + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol, atol = 1.0e-4)) + + def test_df_dh(self): + exact = self.model.df_dh(self.s, self.h, self.T) + numer = utility.batch_differentiate(lambda x: self.model.f(self.s, x, self.T), self.h) + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + + def test_dmix_fn_dh(self): + exact = self.model.dmix_fn_dh(self.s, self.h, self.T) + numer = utility.batch_differentiate(lambda x: self.model.mix_fn(self.s, x, self.T), self.h) + + self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) diff --git a/test/test_models.py b/test/test_models.py index 61e87a0..148a56d 100644 --- a/test/test_models.py +++ b/test/test_models.py @@ -71,7 +71,12 @@ def test_derivs_stress(self): ddv = utility.batch_differentiate( lambda x: use.forward(self.t, x)[0], self.state_stress ) - + + # 3 time steps, 3 state variables + print(stress_rates.shape) + print(dv[1]) + print(ddv[1]) + self.assertTrue(np.allclose(dv, ddv, rtol=1e-4, atol=1e-4)) def test_partial_state(self): @@ -114,6 +119,9 @@ def test_partial_erate(self): lambda y: self.model.forward(t, self.state_strain, y, T)[0], erate ).squeeze(-1) + print(exact) + print(numer) + self.assertTrue(torch.allclose(exact, numer, atol=1e-4, rtol=1e-4)) @@ -330,13 +338,13 @@ def setUp(self): class TestIsoKinRIPlasticity(unittest.TestCase, CommonModel, CommonModelBatchBatch): def setUp(self): self.E = torch.tensor(100000.0) - self.sy = torch.tensor(10.0) + self.sy = torch.tensor(100.0) self.R = torch.tensor(101.0) self.d = torch.tensor(1.3) self.iso = hardening.VoceIsotropicHardeningModel(CP(self.R), CP(self.d)) - self.C = torch.tensor(12000.0) + self.C = torch.tensor(1200.0) self.g = torch.tensor(10.1) self.kin = hardening.FAKinematicHardeningModel(CP(self.C), CP(self.g)) @@ -364,7 +372,7 @@ def setUp(self): ) / 3.0 self.state_stress = ( torch.tensor([[0.05, 30.0, 10.0], [0.07, 10.0, 15.0], [0.08, 50.0, 60.0]]) - ) / 3.0 + ) self.t = self.times[2] self.step = 2 From ed0c0a1b7ffef4c1d27da906e72fb1968318599b Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 11 Jan 2024 20:16:59 -0600 Subject: [PATCH 10/18] Thank god --- pyoptmat/models.py | 8 +++++++- test/test_models.py | 8 -------- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 09f6660..282e17e 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -563,10 +563,16 @@ def RJ(erate): ydot, J, Je, _ = self.model(t, torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), erate, cT) + + # There is an annoying extra term that is the derivative of the history rate with respect to the + # solved for strain rate times the derivative of the strain rate with respect to history + t1 = Je[...,1:].unsqueeze(-1) + t2 = utility.mbmm(1.0/Je[...,:1].unsqueeze(-1), J[...,0,1:].unsqueeze(-2)) + t3 = utility.mbmm(t1, t2) # Corrected jacobian row1 = torch.cat([torch.zeros_like(J[...,0,0]).unsqueeze(-1), -J[..., 0, 1:] / Je[...,0][...,None]], dim = -1).unsqueeze(-2) - rest = torch.cat([torch.zeros_like(J[...,1:,0]).unsqueeze(-1), J[...,1:,1:]], dim = -1) + rest = torch.cat([torch.zeros_like(J[...,1:,0]).unsqueeze(-1), J[...,1:,1:] - t3], dim = -1) jac = torch.cat([row1,rest], dim = -2) return torch.cat([erate.unsqueeze(-1), ydot[...,1:]], dim = -1), jac diff --git a/test/test_models.py b/test/test_models.py index 148a56d..8db3066 100644 --- a/test/test_models.py +++ b/test/test_models.py @@ -72,11 +72,6 @@ def test_derivs_stress(self): lambda x: use.forward(self.t, x)[0], self.state_stress ) - # 3 time steps, 3 state variables - print(stress_rates.shape) - print(dv[1]) - print(ddv[1]) - self.assertTrue(np.allclose(dv, ddv, rtol=1e-4, atol=1e-4)) def test_partial_state(self): @@ -119,9 +114,6 @@ def test_partial_erate(self): lambda y: self.model.forward(t, self.state_strain, y, T)[0], erate ).squeeze(-1) - print(exact) - print(numer) - self.assertTrue(torch.allclose(exact, numer, atol=1e-4, rtol=1e-4)) From ea9f5825d0e135b481202fc6b69b51620d296158 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 11 Jan 2024 21:26:34 -0600 Subject: [PATCH 11/18] Updated tests and blacked --- examples/ode/damping.py | 255 +++++++++----- examples/ode/neuron.py | 324 ++++++++++-------- examples/ode/trajectory.py | 12 +- examples/ode/vanderpohl.py | 44 ++- examples/performance/linear_solve_profile.py | 63 ++-- .../creep-fatigue/deterministic/optimize.py | 6 +- .../creep-fatigue/statistical/infer.py | 4 +- .../tension-temperature/generate_data.py | 108 +++--- .../tension-temperature/infer.py | 145 +++++--- .../tension-temperature/survey-data.py | 6 +- .../tension/deterministic/optimize.py | 8 +- .../tension/statistical/demo_plot.py | 4 +- .../tension/statistical/infer.py | 17 +- examples/structural-material-models/damage.py | 43 ++- .../kocks-mecking-differentiable.py | 63 +++- .../kocks-mecking.py | 16 +- .../rate_independent.py | 25 +- .../rate_independent_approximations.py | 125 ------- .../rate_independent_stability.py | 123 ++++--- .../stress_strain_control.py | 6 +- pyoptmat/chunktime.py | 43 ++- pyoptmat/flowrules.py | 149 +++++--- pyoptmat/models.py | 137 +++++--- pyoptmat/ode.py | 27 +- pyoptmat/optimize.py | 17 +- pyoptmat/solvers.py | 40 ++- test/test_flowrules.py | 87 +++-- test/test_models.py | 8 +- test/test_ode.py | 5 +- 29 files changed, 1130 insertions(+), 780 deletions(-) delete mode 100755 examples/structural-material-models/rate_independent_approximations.py diff --git a/examples/ode/damping.py b/examples/ode/damping.py index 07a6edb..61bcc7f 100755 --- a/examples/ode/damping.py +++ b/examples/ode/damping.py @@ -31,6 +31,7 @@ import matplotlib.pyplot as plt import sys + sys.path.append("../..") from pyoptmat import ode, utility, neuralode @@ -43,6 +44,7 @@ # I'm assuming single precision will be fine for this but for now use doubles torch.set_default_tensor_type(torch.DoubleTensor) + class MassDamperSpring(torch.nn.Module): """ Mass, spring, damper system of arbitrary size @@ -57,7 +59,8 @@ class MassDamperSpring(torch.nn.Module): M (torch.tensor): (n_chain,) vector of masses force (function): gives f(t) """ - def __init__(self, K, C, M, force, C_scale = 1.0e-4, M_scale = 1.0e-5): + + def __init__(self, K, C, M, force, C_scale=1.0e-4, M_scale=1.0e-5): super().__init__() self.K = torch.nn.Parameter(K) @@ -67,7 +70,7 @@ def __init__(self, K, C, M, force, C_scale = 1.0e-4, M_scale = 1.0e-5): self.C_scale = C_scale self.M_scale = M_scale - + self.half_size = K.shape[0] self.size = 2 * self.half_size @@ -84,63 +87,78 @@ def forward(self, t, y): J: (nchunk, nbatch, size, size) tensor of Jacobians """ ydot = torch.zeros_like(y) - J = torch.zeros(y.shape + (self.size,), device = t.device) + J = torch.zeros(y.shape + (self.size,), device=t.device) # Separate out for convenience - u = y[...,:self.half_size] - v = y[...,self.half_size:] - + u = y[..., : self.half_size] + v = y[..., self.half_size :] + # Differences - du = torch.diff(u, dim = -1) - dv = torch.diff(v, dim = -1) + du = torch.diff(u, dim=-1) + dv = torch.diff(v, dim=-1) # Scaled properties M = self.M * self.M_scale C = self.C * self.C_scale - + # Rate # Velocity - ydot[...,:self.half_size] = v - + ydot[..., : self.half_size] = v + # Springs - ydot[...,self.half_size:-1] += self.K[...,:-1] * du / M[...,:-1] - ydot[...,self.half_size+1:] += -self.K[...,:-1] * du / M[...,1:] - ydot[...,-1] += -self.K[...,-1] * u[...,-1] / M[...,-1] - ydot[...,self.half_size] += self.force(t) / M[...,0] + ydot[..., self.half_size : -1] += self.K[..., :-1] * du / M[..., :-1] + ydot[..., self.half_size + 1 :] += -self.K[..., :-1] * du / M[..., 1:] + ydot[..., -1] += -self.K[..., -1] * u[..., -1] / M[..., -1] + ydot[..., self.half_size] += self.force(t) / M[..., 0] # Dampers - ydot[...,self.half_size:-1] += C[...,:-1] * dv / M[...,:-1] - ydot[...,self.half_size+1:] += -C[...,:-1] * dv / M[...,1:] - ydot[...,-1] += -C[...,-1] * v[...,-1] / M[...,-1] + ydot[..., self.half_size : -1] += C[..., :-1] * dv / M[..., :-1] + ydot[..., self.half_size + 1 :] += -C[..., :-1] * dv / M[..., 1:] + ydot[..., -1] += -C[..., -1] * v[..., -1] / M[..., -1] # Derivative - # Velocity - J[...,:self.half_size,self.half_size:] += torch.eye(self.half_size, device = t.device) + # Velocity + J[..., : self.half_size, self.half_size :] += torch.eye( + self.half_size, device=t.device + ) # Springs - J[...,self.half_size:,:self.half_size] += torch.diag_embed(-self.K / M) - J[...,self.half_size+1:,1:self.half_size] += torch.diag_embed(-self.K[...,:-1] / M[...,1:]) - J[...,self.half_size:,:self.half_size] += torch.diag_embed(self.K[...,:-1] / M[...,:-1], offset = 1) - J[...,self.half_size:,:self.half_size] += torch.diag_embed(self.K[...,:-1] / M[...,1:], offset = -1) + J[..., self.half_size :, : self.half_size] += torch.diag_embed(-self.K / M) + J[..., self.half_size + 1 :, 1 : self.half_size] += torch.diag_embed( + -self.K[..., :-1] / M[..., 1:] + ) + J[..., self.half_size :, : self.half_size] += torch.diag_embed( + self.K[..., :-1] / M[..., :-1], offset=1 + ) + J[..., self.half_size :, : self.half_size] += torch.diag_embed( + self.K[..., :-1] / M[..., 1:], offset=-1 + ) # Dampers - J[...,self.half_size:,self.half_size:] += torch.diag_embed(-C / M) - J[...,self.half_size+1:,self.half_size+1:] += torch.diag_embed(-C[...,:-1] / M[...,1:]) - J[...,self.half_size:,self.half_size:] += torch.diag_embed(C[...,:-1] / M[...,:-1], offset = 1) - J[...,self.half_size:,self.half_size:] += torch.diag_embed(C[...,:-1] / M[...,1:], offset = -1) + J[..., self.half_size :, self.half_size :] += torch.diag_embed(-C / M) + J[..., self.half_size + 1 :, self.half_size + 1 :] += torch.diag_embed( + -C[..., :-1] / M[..., 1:] + ) + J[..., self.half_size :, self.half_size :] += torch.diag_embed( + C[..., :-1] / M[..., :-1], offset=1 + ) + J[..., self.half_size :, self.half_size :] += torch.diag_embed( + C[..., :-1] / M[..., 1:], offset=-1 + ) return ydot, J def initial_condition(self, nsamples): """ - I don't want to learn the initial condition (we could) + I don't want to learn the initial condition (we could) so just return deterministic zeros """ - return torch.zeros(nsamples, self.size, device = device) + return torch.zeros(nsamples, self.size, device=device) + def neural_ode_factory(n_hidden, n_layers, n_inter): """Make a neural ODE to try - + Args: n_hidden (int): number of hidden state variables n_layers (int): network depth @@ -150,31 +168,39 @@ def neural_ode_factory(n_hidden, n_layers, n_inter): n_out = n_hidden + 1 n_inter += 1 - mods = [nn.Linear(n_in, n_inter), nn.ReLU() - ] + list(itertools.chain(*[[nn.Linear(n_inter, n_inter), nn.ReLU()] for i in range(n_layers)]) - ) + [nn.Linear(n_inter, n_out)] + mods = ( + [nn.Linear(n_in, n_inter), nn.ReLU()] + + list( + itertools.chain( + *[[nn.Linear(n_inter, n_inter), nn.ReLU()] for i in range(n_layers)] + ) + ) + + [nn.Linear(n_inter, n_out)] + ) return nn.Sequential(*mods) + def random_walk(time, mean, scale, mag): """ Simulate a random walk with velocity of the given mean and scale """ results = torch.zeros_like(time) for i in range(1, time.shape[0]): - dt = time[i] - time[i-1] - v = mag*torch.normal(mean, scale) - results[i] = results[i-1] + v * dt + dt = time[i] - time[i - 1] + v = mag * torch.normal(mean, scale) + results[i] = results[i - 1] + v * dt return results + if __name__ == "__main__": # Basic parameters - n_chain = 5 # Number of spring-dashpot-mass elements - n_samples = 5 # Number of random force samples - n_time = 512+1 # Number of time steps - integration_method = 'backward-euler' - direct_solver = "thomas" # Batched, block, bidiagonal direct solver method + n_chain = 5 # Number of spring-dashpot-mass elements + n_samples = 5 # Number of random force samples + n_time = 512 + 1 # Number of time steps + integration_method = "backward-euler" + direct_solver = "thomas" # Batched, block, bidiagonal direct solver method # Time chunking -- best value may vary on your system n_chunk = 2**7 @@ -183,22 +209,31 @@ def random_walk(time, mean, scale, mag): t_end = 1.0 # Mean and standard deviation for random force "velocity" - force_mean = torch.tensor(0.0, device = device) - force_std = torch.tensor(1.0, device = device) - force_mag = torch.tensor(1.0, device = device) + force_mean = torch.tensor(0.0, device=device) + force_std = torch.tensor(1.0, device=device) + force_mag = torch.tensor(1.0, device=device) # True values of the mass, damping, and stiffness - K = torch.rand(n_chain, device = device) - C = torch.rand(n_chain, device = device) - M = torch.rand(n_chain, device = device) + K = torch.rand(n_chain, device=device) + C = torch.rand(n_chain, device=device) + M = torch.rand(n_chain, device=device) # Time values - time = torch.linspace(0, t_end, n_time, device = device).unsqueeze(-1).expand(n_time, n_samples) + time = ( + torch.linspace(0, t_end, n_time, device=device) + .unsqueeze(-1) + .expand(n_time, n_samples) + ) # 1. Generate the data # Generate the force values - force = random_walk(time, force_mean.unsqueeze(0).expand(n_samples), force_std.unsqueeze(0).expand(n_samples), force_mag.unsqueeze(0).expand(n_samples)) + force = random_walk( + time, + force_mean.unsqueeze(0).expand(n_samples), + force_std.unsqueeze(0).expand(n_samples), + force_mag.unsqueeze(0).expand(n_samples), + ) # Plot them plt.figure() @@ -216,43 +251,56 @@ def random_walk(time, mean, scale, mag): # Generate the data with torch.no_grad(): - y_data = ode.odeint(model, y0, time, method = integration_method, - block_size = n_chunk, - direct_solve_method = direct_solver) + y_data = ode.odeint( + model, + y0, + time, + method=integration_method, + block_size=n_chunk, + direct_solve_method=direct_solver, + ) # The observations will just be the first entry - observable = y_data[...,0] - + observable = y_data[..., 0] + # Plot them plt.figure() plt.plot(time.cpu().numpy(), observable.cpu().numpy()) plt.show() - + # 2. Infer with the actual model - ode_model = MassDamperSpring(torch.rand(n_chain, device = device), - torch.rand(n_chain, device = device), - torch.rand(n_chain, device = device), force_continuous) - + ode_model = MassDamperSpring( + torch.rand(n_chain, device=device), + torch.rand(n_chain, device=device), + torch.rand(n_chain, device=device), + force_continuous, + ) + # Training parameters... niter = 1000 lr = 1.0e-3 - loss = torch.nn.MSELoss(reduction = "sum") + loss = torch.nn.MSELoss(reduction="sum") # Optimization setup optim = torch.optim.Adam(ode_model.parameters(), lr) + def closure(): optim.zero_grad() - pred = ode.odeint_adjoint(ode_model, y0, time, - method = integration_method, - block_size = n_chunk, - direct_solve_method = direct_solver) - obs = pred[...,0] + pred = ode.odeint_adjoint( + ode_model, + y0, + time, + method=integration_method, + block_size=n_chunk, + direct_solve_method=direct_solver, + ) + obs = pred[..., 0] lossv = loss(obs, observable) lossv.backward() return lossv # Optimization loop - t = tqdm.tqdm(range(niter), total = niter, desc = "Loss: ") + t = tqdm.tqdm(range(niter), total=niter, desc="Loss: ") loss_history = [] for i in t: closs = optim.step(closure) @@ -266,44 +314,54 @@ def closure(): # Make the final predictions with torch.no_grad(): - ode_preds = ode.odeint(ode_model, y0, time, method = integration_method, block_size = n_chunk) - ode_obs = ode_preds[...,0] + ode_preds = ode.odeint( + ode_model, y0, time, method=integration_method, block_size=n_chunk + ) + ode_obs = ode_preds[..., 0] # Make a plot plt.figure() for i in range(n_samples): - l, = plt.plot(time.cpu().numpy()[:,i], observable.cpu().numpy()[:,i]) - plt.plot(time.cpu().numpy()[:,i], ode_obs.cpu().numpy()[:,i], ls = '--', color = l.get_color()) + (l,) = plt.plot(time.cpu().numpy()[:, i], observable.cpu().numpy()[:, i]) + plt.plot( + time.cpu().numpy()[:, i], + ode_obs.cpu().numpy()[:, i], + ls="--", + color=l.get_color(), + ) plt.show() # 2. Infer with a neural ODE # Setup the model - n_hidden = n_chain # I don't know, seems reasonable + n_hidden = n_chain # I don't know, seems reasonable n_layers = 3 - n_inter = n_chain # Same thing, seems reasonable - - nn = neural_ode_factory(n_hidden, n_layers, n_inter).to(device) - y0 = torch.zeros(n_samples, n_hidden+1, device = device) + n_inter = n_chain # Same thing, seems reasonable + + nn = neural_ode_factory(n_hidden, n_layers, n_inter).to(device) + y0 = torch.zeros(n_samples, n_hidden + 1, device=device) - nn_model = neuralode.NeuralODE(nn, lambda t: force_continuous(t).unsqueeze(-1)) + nn_model = neuralode.NeuralODE(nn, lambda t: force_continuous(t).unsqueeze(-1)) # Training parameters niter = 1000 lr = 1.0e-3 - loss = torch.nn.MSELoss(reduction = "sum") - + loss = torch.nn.MSELoss(reduction="sum") + # Optimization setup optim = torch.optim.Adam(nn_model.parameters(), lr) + def closure(): optim.zero_grad() - pred = ode.odeint_adjoint(nn_model, y0, time, method = integration_method, block_size = n_chunk) - obs = pred[...,0] + pred = ode.odeint_adjoint( + nn_model, y0, time, method=integration_method, block_size=n_chunk + ) + obs = pred[..., 0] lossv = loss(obs, observable) lossv.backward() return lossv # Optimization loop - t = tqdm.tqdm(range(niter), total = niter, desc = "Loss: ") + t = tqdm.tqdm(range(niter), total=niter, desc="Loss: ") loss_history = [] for i in t: closs = optim.step(closure) @@ -317,20 +375,37 @@ def closure(): # Make the final predictions with torch.no_grad(): - nn_preds = ode.odeint(nn_model, y0, time, method = integration_method, block_size = n_chunk) - nn_obs = nn_preds[...,0] + nn_preds = ode.odeint( + nn_model, y0, time, method=integration_method, block_size=n_chunk + ) + nn_obs = nn_preds[..., 0] # Make a plot plt.figure() for i in range(n_samples): - l, = plt.plot(time.cpu().numpy()[:,i], observable.cpu().numpy()[:,i]) - plt.plot(time.cpu().numpy()[:,i], nn_obs.cpu().numpy()[:,i], ls = '--', color = l.get_color()) + (l,) = plt.plot(time.cpu().numpy()[:, i], observable.cpu().numpy()[:, i]) + plt.plot( + time.cpu().numpy()[:, i], + nn_obs.cpu().numpy()[:, i], + ls="--", + color=l.get_color(), + ) plt.show() - + # Compare all three plt.figure() for i in range(n_samples): - l, = plt.plot(time.cpu().numpy()[:,i], observable.cpu().numpy()[:,i]) - plt.plot(time.cpu().numpy()[:,i], ode_obs.cpu().numpy()[:,i], ls = '--', color = l.get_color()) - plt.plot(time.cpu().numpy()[:,i], nn_obs.cpu().numpy()[:,i], ls = ':', color = l.get_color()) + (l,) = plt.plot(time.cpu().numpy()[:, i], observable.cpu().numpy()[:, i]) + plt.plot( + time.cpu().numpy()[:, i], + ode_obs.cpu().numpy()[:, i], + ls="--", + color=l.get_color(), + ) + plt.plot( + time.cpu().numpy()[:, i], + nn_obs.cpu().numpy()[:, i], + ls=":", + color=l.get_color(), + ) plt.show() diff --git a/examples/ode/neuron.py b/examples/ode/neuron.py index 7f745ef..352e4bc 100755 --- a/examples/ode/neuron.py +++ b/examples/ode/neuron.py @@ -13,7 +13,8 @@ """ import sys -sys.path.append('../..') + +sys.path.append("../..") import torch import torch.nn as nn @@ -33,87 +34,108 @@ dev = "cpu" device = torch.device(dev) -def driving_current(I_max, R, rate, thold, chold, Ncycles, nload, - nhold, neq): - """ - Setup cyclic current histories - Args: - I_max (tuple): (min,max) uniform distribution of max currents - R (tuple): (min,max) uniform distribution of load ratios - rate (tuple): (min, max) uniform distribution of current rates - thold (tuple): (min, max) uniform distribution of hold at max current times - chold (tuple): (min, max) uniform distribution of hold at min current times - Ncycles (int): number of load cycles - nload (int): number of time steps for current ramp - nhold (int): number of time steps for current holds - neq (int): number of independent current histories +def driving_current(I_max, R, rate, thold, chold, Ncycles, nload, nhold, neq): """ - Ntotal = (4*nload + 2*nhold) * Ncycles + 1 - times = torch.empty((Ntotal,neq)) - I = torch.empty((Ntotal,neq)) + Setup cyclic current histories + + Args: + I_max (tuple): (min,max) uniform distribution of max currents + R (tuple): (min,max) uniform distribution of load ratios + rate (tuple): (min, max) uniform distribution of current rates + thold (tuple): (min, max) uniform distribution of hold at max current times + chold (tuple): (min, max) uniform distribution of hold at min current times + Ncycles (int): number of load cycles + nload (int): number of time steps for current ramp + nhold (int): number of time steps for current holds + neq (int): number of independent current histories + """ + Ntotal = (4 * nload + 2 * nhold) * Ncycles + 1 + times = torch.empty((Ntotal, neq)) + I = torch.empty((Ntotal, neq)) for i in range(neq): - cycle = experiments.generate_random_cycle( - I_max, R, rate, thold, chold) + cycle = experiments.generate_random_cycle(I_max, R, rate, thold, chold) ti, Ii, ci = experiments.sample_cycle_normalized_times( - cycle, Ncycles, nload, nhold) - times[:,i] = torch.tensor(ti) - I[:,i] = torch.tensor(Ii) - - return utility.ArbitraryBatchTimeSeriesInterpolator( - times.to(device), I.to(device)), times.to(device), I.to(device) + cycle, Ncycles, nload, nhold + ) + times[:, i] = torch.tensor(ti) + I[:, i] = torch.tensor(Ii) + + return ( + utility.ArbitraryBatchTimeSeriesInterpolator(times.to(device), I.to(device)), + times.to(device), + I.to(device), + ) + class HodgkinHuxleyCoupledNeurons(torch.nn.Module): """ - As described by :cite:`schwemmer2012theory` + As described by :cite:`schwemmer2012theory` - .. math:: + .. math:: - F(y) = \\begin{bmatrix} F_1(y_1) - \\\\ F_2(y_2) - \\\\ \\vdots - \\\\ F_n(y_n) - \\end{bmatrix} + F(y) = \\begin{bmatrix} F_1(y_1) + \\\\ F_2(y_2) + \\\\ \\vdots + \\\\ F_n(y_n) + \\end{bmatrix} - with + with - .. math:: + .. math:: - y = \\begin{bmatrix} y_1 - \\\ y_2 - \\\ \\vdots - \\\ y_n - \\end{bmatrix} + y = \\begin{bmatrix} y_1 + \\\ y_2 + \\\ \\vdots + \\\ y_n + \\end{bmatrix} - and + and - .. math:: + .. math:: - F_i(y_i) = \\begin{bmatrix} V_i - \\\\ m_i - \\\\ h_i - \\\\ n_i - \\end{bmatrix} + F_i(y_i) = \\begin{bmatrix} V_i + \\\\ m_i + \\\\ h_i + \\\\ n_i + \\end{bmatrix} - and + and - .. math:: + .. math:: - F_i(y_i) = \\begin{bmatrix} \\frac{1}{C_i}(-g_{Na,i}m_i^3h_i(V_i-E_{Na,i})-g_{K,i}n_i^4(V_i-E_{k,i})-g_{L,i}(V_i-E_{L,i})+I(t) + \\Delta V_{ij}) - \\\\ \\frac{m_{\\infty,i}-m_i}{\\tau_{m,i}} - \\\\ \\frac{h_{\\infty,i}-h_i}{\\tau_{h,i}} - \\\\ \\frac{n_{\\infty,i}-n_i}{\\tau_{n,i}} - \\end{bmatrix} + F_i(y_i) = \\begin{bmatrix} \\frac{1}{C_i}(-g_{Na,i}m_i^3h_i(V_i-E_{Na,i})-g_{K,i}n_i^4(V_i-E_{k,i})-g_{L,i}(V_i-E_{L,i})+I(t) + \\Delta V_{ij}) + \\\\ \\frac{m_{\\infty,i}-m_i}{\\tau_{m,i}} + \\\\ \\frac{h_{\\infty,i}-h_i}{\\tau_{h,i}} + \\\\ \\frac{n_{\\infty,i}-n_i}{\\tau_{n,i}} + \\end{bmatrix} - with + with - .. math:: + .. math:: - \\Delta V_{ij} = \\sum_{j=1}^{n_{neurons}}\\frac{g_{C,i}(V_i-V_j)}{C_i} + \\Delta V_{ij} = \\sum_{j=1}^{n_{neurons}}\\frac{g_{C,i}(V_i-V_j)}{C_i} """ - def __init__(self, C, g_Na, E_Na, g_K, E_K, g_L, E_L, m_inf, - tau_m, h_inf, tau_h, n_inf, tau_n, g_C, I): + + def __init__( + self, + C, + g_Na, + E_Na, + g_K, + E_K, + g_L, + E_L, + m_inf, + tau_m, + h_inf, + tau_h, + n_inf, + tau_n, + g_C, + I, + ): super().__init__() self.C = nn.Parameter(C) self.g_Na = nn.Parameter(g_Na) @@ -137,79 +159,106 @@ def __init__(self, C, g_Na, E_Na, g_K, E_K, g_L, E_L, m_inf, def forward(self, t, y): """ - Evaluate the system of ODEs at a particular state + Evaluate the system of ODEs at a particular state - Args: - t (torch.tensor): current time - y (torch.tensor): current state + Args: + t (torch.tensor): current time + y (torch.tensor): current state - Returns: - y_dot (torch.tensor): current ODEs - J (torch.tensor): current jacobian + Returns: + y_dot (torch.tensor): current ODEs + J (torch.tensor): current jacobian """ Ic = self.I(t) - V = y[...,0::4].clone() - m = y[...,1::4].clone() - h = y[...,2::4].clone() - n = y[...,3::4].clone() + V = y[..., 0::4].clone() + m = y[..., 1::4].clone() + h = y[..., 2::4].clone() + n = y[..., 3::4].clone() ydot = torch.zeros_like(y) - + # V - ydot[...,0::4] = 1.0 / self.C[None,...] * ( - -self.g_Na[None,...] * m**3.0 * h * (V - self.E_Na[None,...]) - -self.g_K[None,...] * n**4.0 * (V - self.E_K[None,...]) - -self.g_L[None,...] * (V - self.E_L[None,...]) + Ic[...,None]) + ydot[..., 0::4] = ( + 1.0 + / self.C[None, ...] + * ( + -self.g_Na[None, ...] * m**3.0 * h * (V - self.E_Na[None, ...]) + - self.g_K[None, ...] * n**4.0 * (V - self.E_K[None, ...]) + - self.g_L[None, ...] * (V - self.E_L[None, ...]) + + Ic[..., None] + ) + ) # Coupling term - dV = torch.sum(self.g_C[None,...]*(V[...,:,None] - V[...,None,:]) / self.C[None,...], dim = -1) - ydot[...,0::4] += dV + dV = torch.sum( + self.g_C[None, ...] + * (V[..., :, None] - V[..., None, :]) + / self.C[None, ...], + dim=-1, + ) + ydot[..., 0::4] += dV # m - ydot[...,1::4] = (self.m_inf - m) / self.tau_m + ydot[..., 1::4] = (self.m_inf - m) / self.tau_m - # h - ydot[...,2::4] = (self.h_inf - h) / self.tau_h + # h + ydot[..., 2::4] = (self.h_inf - h) / self.tau_h # n - ydot[...,3::4] = (self.n_inf - n) / self.tau_n + ydot[..., 3::4] = (self.n_inf - n) / self.tau_n + + J = torch.zeros(y.shape + y.shape[-1:], device=ydot.device) - J = torch.zeros(y.shape + y.shape[-1:], device = ydot.device) - # V, V - J[...,0::4,0::4] = torch.diag_embed(1.0 / self.C[None,...] * ( - -self.g_L[None,...] - -self.g_Na[None,...]*h*m**3.0 - -self.g_K[None,...]*n**4.0)) + J[..., 0::4, 0::4] = torch.diag_embed( + 1.0 + / self.C[None, ...] + * ( + -self.g_L[None, ...] + - self.g_Na[None, ...] * h * m**3.0 + - self.g_K[None, ...] * n**4.0 + ) + ) # Coupling term - J[...,0::4,0::4] -= self.g_C[None,...] / self.C[None,...] - - J[...,0::4,0::4] += torch.eye(self.n_neurons, device = ydot.device).expand( - *ydot.shape[:-1], -1, -1) * torch.sum(self.g_C / self.C) - + J[..., 0::4, 0::4] -= self.g_C[None, ...] / self.C[None, ...] + + J[..., 0::4, 0::4] += torch.eye(self.n_neurons, device=ydot.device).expand( + *ydot.shape[:-1], -1, -1 + ) * torch.sum(self.g_C / self.C) + # V, m - J[...,0::4,1::4] = torch.diag_embed(-1.0 / self.C[None,...] * ( - 3.0 * self.g_Na[None,...] * h * m**2.0 * (-self.E_Na[None,...] + V))) + J[..., 0::4, 1::4] = torch.diag_embed( + -1.0 + / self.C[None, ...] + * (3.0 * self.g_Na[None, ...] * h * m**2.0 * (-self.E_Na[None, ...] + V)) + ) # V, h - J[...,0::4,2::4] = torch.diag_embed(-1.0 / self.C[None,...] * ( - self.g_Na[None,...] * m**3.0 * (-self.E_Na[None,...] + V))) + J[..., 0::4, 2::4] = torch.diag_embed( + -1.0 + / self.C[None, ...] + * (self.g_Na[None, ...] * m**3.0 * (-self.E_Na[None, ...] + V)) + ) # V, n - J[...,0::4,3::4] = torch.diag_embed(-1.0 / self.C[None,...] * ( - 4.0 * self.g_K[None,...] * n**3.0 * (-self.E_K[None,...] + V))) + J[..., 0::4, 3::4] = torch.diag_embed( + -1.0 + / self.C[None, ...] + * (4.0 * self.g_K[None, ...] * n**3.0 * (-self.E_K[None, ...] + V)) + ) # m, m - J[...,1::4,1::4] = torch.diag(-1.0 / self.tau_m) + J[..., 1::4, 1::4] = torch.diag(-1.0 / self.tau_m) - # h, h - J[...,2::4,2::4] = torch.diag(-1.0 / self.tau_h) + # h, h + J[..., 2::4, 2::4] = torch.diag(-1.0 / self.tau_h) # n, n - J[...,3::4,3::4] = torch.diag(-1.0 / self.tau_n) - + J[..., 3::4, 3::4] = torch.diag(-1.0 / self.tau_n) + return ydot, J + if __name__ == "__main__": # Batch size nbatch = 100 @@ -221,56 +270,55 @@ def forward(self, t, y): N = 5 # Number of time steps per cycle n = 40 - + # Setup the driving current - current, times, discrete_current = driving_current([0.1,1], - [-1,0.9], - [0.5,1.5], - [0,1], - [0,1], - N, - n, - n, - nbatch) - - plt.plot(times.cpu().numpy()[:,::4], discrete_current.cpu().numpy()[:,::4]) + current, times, discrete_current = driving_current( + [0.1, 1], [-1, 0.9], [0.5, 1.5], [0, 1], [0, 1], N, n, n, nbatch + ) + + plt.plot(times.cpu().numpy()[:, ::4], discrete_current.cpu().numpy()[:, ::4]) plt.xlabel("Time") plt.ylabel("Current") plt.show() # Setup model and initial conditions nsteps = times.shape[0] - + model = HodgkinHuxleyCoupledNeurons( - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device), - torch.rand((neq,), device = device) * 5.0, - torch.rand((neq,), device = device), - torch.rand((neq,), device = device) * 15.0, - torch.rand((neq,), device = device), - torch.rand((neq,), device = device) * 10.0, - torch.rand((neq,), device = device) * 0.01, - current) - - y0 = torch.rand((nbatch,model.n_equations), device = device) - + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device), + torch.rand((neq,), device=device) * 5.0, + torch.rand((neq,), device=device), + torch.rand((neq,), device=device) * 15.0, + torch.rand((neq,), device=device), + torch.rand((neq,), device=device) * 10.0, + torch.rand((neq,), device=device) * 0.01, + current, + ) + + y0 = torch.rand((nbatch, model.n_equations), device=device) # Include a backward pass to give a better example of timing t1 = time.time() - res_imp = ode.odeint_adjoint(model, y0, times, block_size = time_block) + res_imp = ode.odeint_adjoint(model, y0, times, block_size=time_block) loss = torch.norm(res_imp) loss.backward() etime = time.time() - t1 - - print("%i batch size, %i blocked time steps, %i neurons, %i time steps: %f s" % (nbatch, time_block, neq, nsteps, etime)) - plt.plot(times[:,0].detach().cpu().numpy(), res_imp[:,0,0::4].detach().cpu().numpy()) + print( + "%i batch size, %i blocked time steps, %i neurons, %i time steps: %f s" + % (nbatch, time_block, neq, nsteps, etime) + ) + + plt.plot( + times[:, 0].detach().cpu().numpy(), res_imp[:, 0, 0::4].detach().cpu().numpy() + ) plt.xlabel("Time") plt.ylabel("Voltage") plt.show() diff --git a/examples/ode/trajectory.py b/examples/ode/trajectory.py index fbf1e50..26a55b5 100755 --- a/examples/ode/trajectory.py +++ b/examples/ode/trajectory.py @@ -63,6 +63,7 @@ # Prior for the noise eps_prior = 0.1 # Just measure variance in data... + def model_act(times): """ times: ntime x nbatch @@ -77,7 +78,7 @@ def model_act(times): torch.stack( ( v * torch.cos(a) * times, - v * torch.sin(a) * times - 0.5 * g * times ** 2.0, + v * torch.sin(a) * times - 0.5 * g * times**2.0, ) ).T, eps_act, @@ -88,7 +89,7 @@ def model_act(times): class Integrator(pyro.nn.PyroModule): - def __init__(self, eqn, y0, extra_params=[], block_size = 1): + def __init__(self, eqn, y0, extra_params=[], block_size=1): super().__init__() self.eqn = eqn self.y0 = y0 @@ -100,7 +101,7 @@ def forward(self, times): self.eqn, self.y0, times, - block_size = self.block_size, + block_size=self.block_size, extra_params=self.extra_params, ) @@ -274,8 +275,9 @@ def gen_extra(self): pyro.clear_param_store() def maker(v, a, **kwargs): - return Integrator(ODE(v, a), torch.zeros(nsamples, 2), - block_size = time_block, **kwargs) + return Integrator( + ODE(v, a), torch.zeros(nsamples, 2), block_size=time_block, **kwargs + ) # Setup the model model = Model( diff --git a/examples/ode/vanderpohl.py b/examples/ode/vanderpohl.py index f42fe24..94327b1 100755 --- a/examples/ode/vanderpohl.py +++ b/examples/ode/vanderpohl.py @@ -81,11 +81,21 @@ def forward(self, t, y): model = VanderPolODE(mu) - res_exp = ode.odeint(model, y0, times, method="forward-euler", - block_size = time_chunk, guess_type = "previous") + res_exp = ode.odeint( + model, + y0, + times, + method="forward-euler", + block_size=time_chunk, + guess_type="previous", + ) res_imp = ode.odeint( - model, y0, times, method="backward-euler", - block_size = time_chunk, guess_type = "previous" + model, + y0, + times, + method="backward-euler", + block_size=time_chunk, + guess_type="previous", ) plt.figure() @@ -103,11 +113,21 @@ def forward(self, t, y): model = VanderPolODE(mu) - res_exp = ode.odeint(model, y0, times, method="forward-euler", - block_size = time_chunk, guess_type = "previous") + res_exp = ode.odeint( + model, + y0, + times, + method="forward-euler", + block_size=time_chunk, + guess_type="previous", + ) res_imp = ode.odeint( - model, y0, times, method="backward-euler", - block_size = time_chunk, guess_type = "previous" + model, + y0, + times, + method="backward-euler", + block_size=time_chunk, + guess_type="previous", ) plt.figure() @@ -126,8 +146,12 @@ def forward(self, t, y): model = VanderPolODE(mu) res_imp = ode.odeint( - model, y0, times, method="backward-euler", - block_size = time_chunk, guess_type = "previous" + model, + y0, + times, + method="backward-euler", + block_size=time_chunk, + guess_type="previous", ) plt.figure() diff --git a/examples/performance/linear_solve_profile.py b/examples/performance/linear_solve_profile.py index acfabfd..f47088a 100755 --- a/examples/performance/linear_solve_profile.py +++ b/examples/performance/linear_solve_profile.py @@ -21,7 +21,8 @@ dev = "cpu" device = torch.device(dev) -def run_time(operator, D, L, v, repeat = 1): + +def run_time(operator, D, L, v, repeat=1): times = [] for i in range(repeat): t1 = time.time() @@ -31,6 +32,7 @@ def run_time(operator, D, L, v, repeat = 1): return np.mean(times) + if __name__ == "__main__": # Number of repeated trials to average over avg = 3 @@ -42,41 +44,54 @@ def run_time(operator, D, L, v, repeat = 1): max_blk = 1000 # Number of samples in range(1,max_blk) to sample num_samples = 10 - - nblks = sorted(random.sample(list(range(1,max_blk)), num_samples)) - - methods = [chunktime.BidiagonalThomasFactorization, chunktime.BidiagonalPCRFactorization, - lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size = 8), - lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size = 16), - lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size = 32), - lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size = 64), - lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size = 128)] - method_names = ["Thomas", "PCR", "Hybrid, n = 8", "Hybrid, n = 16", "Hybrid, n = 32", - "Hybrid, n = 64", "Hybrid, n = 128"] + nblks = sorted(random.sample(list(range(1, max_blk)), num_samples)) + + methods = [ + chunktime.BidiagonalThomasFactorization, + chunktime.BidiagonalPCRFactorization, + lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size=8), + lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size=16), + lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size=32), + lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size=64), + lambda A, B: chunktime.BidiagonalHybridFactorization(A, B, min_size=128), + ] + + method_names = [ + "Thomas", + "PCR", + "Hybrid, n = 8", + "Hybrid, n = 16", + "Hybrid, n = 32", + "Hybrid, n = 64", + "Hybrid, n = 128", + ] nmethods = len(methods) ncase = len(nblks) times = np.zeros((nmethods, ncase)) - - # Do this once to warm up the GPU, it seems to matter - run_time(methods[0], torch.rand(3, nbat, nsize, nsize, device = device), - torch.rand(2, nbat, nsize, nsize, device = device) / 10.0, - torch.rand(nbat, 3 * nsize, device = device)) - for i,nblk in enumerate(nblks): + # Do this once to warm up the GPU, it seems to matter + run_time( + methods[0], + torch.rand(3, nbat, nsize, nsize, device=device), + torch.rand(2, nbat, nsize, nsize, device=device) / 10.0, + torch.rand(nbat, 3 * nsize, device=device), + ) + + for i, nblk in enumerate(nblks): print(nblk) - D = torch.rand(nblk, nbat, nsize, nsize, device = device) - L = torch.rand(nblk - 1, nbat, nsize, nsize, device = device) / 10.0 + D = torch.rand(nblk, nbat, nsize, nsize, device=device) + L = torch.rand(nblk - 1, nbat, nsize, nsize, device=device) / 10.0 - v = torch.rand(nbat, nblk * nsize, device = device) + v = torch.rand(nbat, nblk * nsize, device=device) for j, method in enumerate(methods): - times[j,i] = run_time(method, D, L, v, repeat = avg) + times[j, i] = run_time(method, D, L, v, repeat=avg) - data = pd.DataFrame(data = times.T, index = nblks, columns = method_names) + data = pd.DataFrame(data=times.T, index=nblks, columns=method_names) data.avg = avg data.nsize = nsize data.nbat = nbat - + data.to_csv(f"{nbat}_{nsize}.csv") diff --git a/examples/structural-inference/creep-fatigue/deterministic/optimize.py b/examples/structural-inference/creep-fatigue/deterministic/optimize.py index 3512903..e6aa158 100755 --- a/examples/structural-inference/creep-fatigue/deterministic/optimize.py +++ b/examples/structural-inference/creep-fatigue/deterministic/optimize.py @@ -76,7 +76,11 @@ def make(n, eta, s0, R, d, C, g, **kwargs): print("") # 3) Create the actual model - model = optimize.DeterministicModel(lambda *args, **kwargs: make(*args, block_size = time_chunk_size, **kwargs), names, ics) + model = optimize.DeterministicModel( + lambda *args, **kwargs: make(*args, block_size=time_chunk_size, **kwargs), + names, + ics, + ) # 4) Setup the optimizer niter = 200 diff --git a/examples/structural-inference/creep-fatigue/statistical/infer.py b/examples/structural-inference/creep-fatigue/statistical/infer.py index a7f28dd..df4f2b3 100755 --- a/examples/structural-inference/creep-fatigue/statistical/infer.py +++ b/examples/structural-inference/creep-fatigue/statistical/infer.py @@ -87,14 +87,14 @@ def make(n, eta, s0, R, d, C, g, **kwargs): # 3) Create the actual model model = optimize.HierarchicalStatisticalModel( - lambda *args, **kwargs: make(*args, block_size = time_chunk_size, **kwargs), + lambda *args, **kwargs: make(*args, block_size=time_chunk_size, **kwargs), names, loc_loc_priors, loc_scale_priors, scale_scale_priors, eps, include_noise=False, - use_cached_guess= True + use_cached_guess=True, ).to(device) # 4) Get the guide diff --git a/examples/structural-inference/tension-temperature/generate_data.py b/examples/structural-inference/tension-temperature/generate_data.py index 1e83819..407da21 100755 --- a/examples/structural-inference/tension-temperature/generate_data.py +++ b/examples/structural-inference/tension-temperature/generate_data.py @@ -5,7 +5,8 @@ """ import sys -sys.path.append('../../..') + +sys.path.append("../../..") import numpy as np import torch @@ -27,8 +28,12 @@ # True model parameters # Elastic constants -E_poly = torch.tensor([-5.76107011e-05, 7.52478382e-02, -9.98116448e+01, 2.19193150e+05], device = device) -mu_poly = torch.tensor([-2.19888172e-05, 2.87205489e-02, -3.80960476e+01, 8.36615078e+04], device = device) +E_poly = torch.tensor( + [-5.76107011e-05, 7.52478382e-02, -9.98116448e01, 2.19193150e05], device=device +) +mu_poly = torch.tensor( + [-2.19888172e-05, 2.87205489e-02, -3.80960476e01, 8.36615078e04], device=device +) # Rate sensitivity A = -3.35 @@ -42,15 +47,16 @@ Q = 800.0 # Various physical/reference constants -eps0 = torch.tensor(1.0e6, device = device) -b = torch.tensor(2.019e-7, device = device) -k = torch.tensor(1.38064e-20, device = device) +eps0 = torch.tensor(1.0e6, device=device) +b = torch.tensor(2.019e-7, device=device) +k = torch.tensor(1.38064e-20, device=device) # Setup constants E = temperature.PolynomialScaling(E_poly) mu = temperature.PolynomialScaling(mu_poly) -def true_model(A, B, H, tau0, Q, device = torch.device("cpu"), **kwargs): + +def true_model(A, B, H, tau0, Q, device=torch.device("cpu"), **kwargs): """True analytic model Args: @@ -65,26 +71,21 @@ def true_model(A, B, H, tau0, Q, device = torch.device("cpu"), **kwargs): theta = temperature.ShearModulusScaling(H, mu) tau = temperature.InverseArrheniusScaling(tau0, Q) - isotropic = hardening.Theta0VoceIsotropicHardeningModel( - tau, - theta - ) + isotropic = hardening.Theta0VoceIsotropicHardeningModel(tau, theta) kinematic = hardening.NoKinematicHardeningModel() flowrule = flowrules.IsoKinViscoplasticity( - n, - eta, - temperature.ConstantParameter(torch.tensor(0.0, device = device)), - isotropic, - kinematic - ) - model = models.InelasticModel( - E, - flowrule + n, + eta, + temperature.ConstantParameter(torch.tensor(0.0, device=device)), + isotropic, + kinematic, ) + model = models.InelasticModel(E, flowrule) return models.ModelIntegrator(model, **kwargs).to(device) + if __name__ == "__main__": # Amount of variability sf = 0.075 @@ -101,7 +102,7 @@ def true_model(A, B, H, tau0, Q, device = torch.device("cpu"), **kwargs): nrate = 5 ntemp = 50 strain_rates = np.logspace(-5, -2, nrate) - temperatures = np.linspace(500+273.15,800+273.15, ntemp) + temperatures = np.linspace(500 + 273.15, 800 + 273.15, ntemp) # Number of steps nsteps = 75 @@ -111,60 +112,75 @@ def true_model(A, B, H, tau0, Q, device = torch.device("cpu"), **kwargs): SR = SR.flatten() T = T.flatten() elimits = np.ones_like(T) * elimit - + # Make the actual input data - times, strains, temps, cycles = experiments.make_tension_tests(SR, T, elimits, nsteps) - data = torch.stack((times, temps,strains)) - control = torch.zeros((len(SR),), dtype = int) + times, strains, temps, cycles = experiments.make_tension_tests( + SR, T, elimits, nsteps + ) + data = torch.stack((times, temps, strains)) + control = torch.zeros((len(SR),), dtype=int) types = torch.zeros_like(control) # Setup the parameters - locs = torch.tensor([A, B, H, tau0, Q], device = device) + locs = torch.tensor([A, B, H, tau0, Q], device=device) scales = sf * torch.abs(locs) - noise = torch.tensor(noise, device = device) + noise = torch.tensor(noise, device=device) # Setup the statistical model names = ["A", "B", "H", "tau0", "Q"] model = optimize.StatisticalModel( - lambda *args, **kwargs: true_model(*args, block_size = time_chunk_size, **kwargs), - names, - locs.to(device), scales.to(device), noise.to(device)) - + lambda *args, **kwargs: true_model(*args, block_size=time_chunk_size, **kwargs), + names, + locs.to(device), + scales.to(device), + noise.to(device), + ) + # Data storage... - full_stress = torch.zeros_like(strains.unsqueeze(0).repeat_interleave(nsamples, dim = 0)) - + full_stress = torch.zeros_like( + strains.unsqueeze(0).repeat_interleave(nsamples, dim=0) + ) + # Run forward samples for i in tqdm.trange(nsamples): with torch.no_grad(): - full_stress[i,...] = model(data.to(device), cycles.to(device), types.to(device), control.to(device)) - + full_stress[i, ...] = model( + data.to(device), cycles.to(device), types.to(device), control.to(device) + ) + # Expand the input data - full_times = times.unsqueeze(0).repeat_interleave(nsamples, dim = 0) - full_strains = strains.unsqueeze(0).repeat_interleave(nsamples, dim = 0) - full_temps = temps.unsqueeze(0).repeat_interleave(nsamples, dim = 0) - full_cycles = cycles.unsqueeze(0).repeat_interleave(nsamples, dim = 0) - full_control = control.unsqueeze(0).repeat_interleave(nsamples, dim = 0) - full_types = types.unsqueeze(0).repeat_interleave(nsamples, dim = 0) + full_times = times.unsqueeze(0).repeat_interleave(nsamples, dim=0) + full_strains = strains.unsqueeze(0).repeat_interleave(nsamples, dim=0) + full_temps = temps.unsqueeze(0).repeat_interleave(nsamples, dim=0) + full_cycles = cycles.unsqueeze(0).repeat_interleave(nsamples, dim=0) + full_control = control.unsqueeze(0).repeat_interleave(nsamples, dim=0) + full_types = types.unsqueeze(0).repeat_interleave(nsamples, dim=0) string_control = np.array(["strain"] * full_control.flatten().shape[0]) string_types = np.array(["tensile"] * full_types.flatten().shape[0]) ds = xr.Dataset( { - "time": (["ntime", "nexp"], full_times.transpose(0,1).flatten(-2,-1).cpu().numpy()), + "time": ( + ["ntime", "nexp"], + full_times.transpose(0, 1).flatten(-2, -1).cpu().numpy(), + ), "strain": ( ["ntime", "nexp"], - full_strains.transpose(0,1).flatten(-2,-1).cpu().numpy(), + full_strains.transpose(0, 1).flatten(-2, -1).cpu().numpy(), ), "stress": ( ["ntime", "nexp"], - full_stress.transpose(0,1).flatten(-2,-1).cpu().numpy(), + full_stress.transpose(0, 1).flatten(-2, -1).cpu().numpy(), ), "temperature": ( ["ntime", "nexp"], - full_temps.transpose(0,1).flatten(-2,-1).cpu().numpy(), + full_temps.transpose(0, 1).flatten(-2, -1).cpu().numpy(), + ), + "cycle": ( + ["ntime", "nexp"], + full_cycles.transpose(0, 1).flatten(-2, -1).cpu().numpy(), ), - "cycle": (["ntime", "nexp"], full_cycles.transpose(0,1).flatten(-2,-1).cpu().numpy()), "type": (["nexp"], string_types), "control": (["nexp"], string_control), }, diff --git a/examples/structural-inference/tension-temperature/infer.py b/examples/structural-inference/tension-temperature/infer.py index 3984e39..14a48d8 100755 --- a/examples/structural-inference/tension-temperature/infer.py +++ b/examples/structural-inference/tension-temperature/infer.py @@ -1,7 +1,8 @@ #!/usr/bin/env python3 import sys -sys.path.append('../../..') + +sys.path.append("../../..") import xarray as xr import torch @@ -10,7 +11,15 @@ import tqdm -from pyoptmat import optimize, experiments, models, flowrules, hardening, temperature, scaling +from pyoptmat import ( + optimize, + experiments, + models, + flowrules, + hardening, + temperature, + scaling, +) # Use doubles torch.set_default_tensor_type(torch.DoubleTensor) @@ -23,7 +32,9 @@ device = torch.device(dev) # Elastic constants, don't infer -E_poly = torch.tensor([-5.76107011e-05, 7.52478382e-02, -9.98116448e+01, 2.19193150e+05], device = device) +E_poly = torch.tensor( + [-5.76107011e-05, 7.52478382e-02, -9.98116448e01, 2.19193150e05], device=device +) E = temperature.PolynomialScaling(E_poly) if __name__ == "__main__": @@ -32,82 +43,111 @@ # Load in the data input_data = xr.open_dataset("data.nc") - data, results, cycles, types, control = experiments.load_results( - input_data) + data, results, cycles, types, control = experiments.load_results(input_data) # Figure out our temperature control points ncontrol = 5 - Tcontrol = torch.linspace(torch.min(data[1]), torch.max(data[1]), - ncontrol, device = device) - + Tcontrol = torch.linspace( + torch.min(data[1]), torch.max(data[1]), ncontrol, device=device + ) + # Model making function - def make(n_vals, eta_vals, theta_vals, tau_vals, scale_functions = [lambda x: x] * 4, **kwargs): - n = temperature.PiecewiseScaling(Tcontrol, n_vals, values_scale_fn = scale_functions[0]) - eta = temperature.PiecewiseScaling(Tcontrol, eta_vals, values_scale_fn = scale_functions[1]) - theta = temperature.PiecewiseScaling(Tcontrol, theta_vals, values_scale_fn = scale_functions[2]) - tau = temperature.PiecewiseScaling(Tcontrol, tau_vals, values_scale_fn=scale_functions[3]) + def make( + n_vals, + eta_vals, + theta_vals, + tau_vals, + scale_functions=[lambda x: x] * 4, + **kwargs + ): + n = temperature.PiecewiseScaling( + Tcontrol, n_vals, values_scale_fn=scale_functions[0] + ) + eta = temperature.PiecewiseScaling( + Tcontrol, eta_vals, values_scale_fn=scale_functions[1] + ) + theta = temperature.PiecewiseScaling( + Tcontrol, theta_vals, values_scale_fn=scale_functions[2] + ) + tau = temperature.PiecewiseScaling( + Tcontrol, tau_vals, values_scale_fn=scale_functions[3] + ) - isotropic = hardening.Theta0VoceIsotropicHardeningModel( - tau, - theta - ) + isotropic = hardening.Theta0VoceIsotropicHardeningModel(tau, theta) kinematic = hardening.NoKinematicHardeningModel() flowrule = flowrules.IsoKinViscoplasticity( - n, - eta, - temperature.ConstantParameter(torch.tensor(0.0, device = device)), - isotropic, - kinematic - ) - model = models.InelasticModel( - E, - flowrule + n, + eta, + temperature.ConstantParameter(torch.tensor(0.0, device=device)), + isotropic, + kinematic, ) + model = models.InelasticModel(E, flowrule) return models.ModelIntegrator(model, **kwargs).to(device) - + # Setup initial guesses names = ["n_vals", "eta_vals", "theta_vals", "tau_vals"] scale_functions = [ - scaling.BoundedScalingFunction(torch.tensor(1.1, device = device), torch.tensor(15.0, device = device)), - scaling.BoundedScalingFunction(torch.tensor(20.0, device = device), torch.tensor(1000.0, device = device)), - scaling.BoundedScalingFunction(torch.tensor(200.0, device = device), torch.tensor(20000.0, device = device)), - scaling.BoundedScalingFunction(torch.tensor(10.0, device = device), torch.tensor(1000.0, device = device)) - ] + scaling.BoundedScalingFunction( + torch.tensor(1.1, device=device), torch.tensor(15.0, device=device) + ), + scaling.BoundedScalingFunction( + torch.tensor(20.0, device=device), torch.tensor(1000.0, device=device) + ), + scaling.BoundedScalingFunction( + torch.tensor(200.0, device=device), torch.tensor(20000.0, device=device) + ), + scaling.BoundedScalingFunction( + torch.tensor(10.0, device=device), torch.tensor(1000.0, device=device) + ), + ] loc_loc_priors = [ - scale_functions[0].unscale(torch.ones(ncontrol, device = device) * 5.0), - scale_functions[1].unscale(torch.ones(ncontrol, device = device) * 500.0), - scale_functions[2].unscale(torch.ones(ncontrol, device = device) * 3000.0), - scale_functions[3].unscale(torch.ones(ncontrol, device = device) * 500.0) - ] + scale_functions[0].unscale(torch.ones(ncontrol, device=device) * 5.0), + scale_functions[1].unscale(torch.ones(ncontrol, device=device) * 500.0), + scale_functions[2].unscale(torch.ones(ncontrol, device=device) * 3000.0), + scale_functions[3].unscale(torch.ones(ncontrol, device=device) * 500.0), + ] loc_scale_priors = [0.1 * l for l in loc_loc_priors] scale_scale_priors = [0.1 * l for l in loc_loc_priors] - eps_prior = torch.tensor(5.0, device = device) + eps_prior = torch.tensor(5.0, device=device) print("Prior values") - for n,l,s1, s2 in zip(names, loc_loc_priors, loc_scale_priors, scale_scale_priors): + for n, l, s1, s2 in zip( + names, loc_loc_priors, loc_scale_priors, scale_scale_priors + ): print(n) print("\tloc: %s" % str(l.cpu())) print("\tloc scale: %s" % str(s1.cpu())) print("\tscale scale: %s" % str(s2.cpu())) # Do a quick check on the reasonableness of the priors - test = optimize.DeterministicModel(lambda *args, **kwargs: make(*args, - scale_functions=scale_functions, block_size = time_chunk_size, **kwargs), - names, loc_loc_priors) - + test = optimize.DeterministicModel( + lambda *args, **kwargs: make( + *args, scale_functions=scale_functions, block_size=time_chunk_size, **kwargs + ), + names, + loc_loc_priors, + ) + with torch.no_grad(): - test_results = test(data.to(device), cycles.to(device), types.to(device), - control.to(device)) + test_results = test( + data.to(device), cycles.to(device), types.to(device), control.to(device) + ) plt.plot(data[-1].cpu().numpy(), test_results.cpu().numpy()) plt.show() # Create the actual model model = optimize.HierarchicalStatisticalModel( - lambda *args, **kwargs: make(*args, scale_functions = scale_functions, - block_size = time_chunk_size, **kwargs), - names, loc_loc_priors, loc_scale_priors, scale_scale_priors, eps_prior + lambda *args, **kwargs: make( + *args, scale_functions=scale_functions, block_size=time_chunk_size, **kwargs + ), + names, + loc_loc_priors, + loc_scale_priors, + scale_scale_priors, + eps_prior, ).to(device) # Get the guide @@ -129,7 +169,12 @@ def make(n_vals, eta_vals, theta_vals, tau_vals, scale_functions = [lambda x: x] t = tqdm.tqdm(range(niter), total=niter, desc="Loss: ") loss_hist = [] for i in t: - loss = svi.step(data.to(device), cycles.to(device), types.to(device), - control.to(device), results.to(device)) + loss = svi.step( + data.to(device), + cycles.to(device), + types.to(device), + control.to(device), + results.to(device), + ) loss_hist.append(loss) t.set_description("Loss %3.2e" % loss) diff --git a/examples/structural-inference/tension-temperature/survey-data.py b/examples/structural-inference/tension-temperature/survey-data.py index ed54a43..8418a95 100755 --- a/examples/structural-inference/tension-temperature/survey-data.py +++ b/examples/structural-inference/tension-temperature/survey-data.py @@ -5,7 +5,8 @@ """ import sys -sys.path.append('../../..') + +sys.path.append("../../..") import numpy as np import torch @@ -19,8 +20,7 @@ if __name__ == "__main__": input_data = xr.open_dataset("data.nc") - data, results, cycles, types, control = experiments.load_results( - input_data) + data, results, cycles, types, control = experiments.load_results(input_data) plt.plot(data[-1].numpy(), results.numpy()) plt.xlabel("Strain (mm/mm)") diff --git a/examples/structural-inference/tension/deterministic/optimize.py b/examples/structural-inference/tension/deterministic/optimize.py index 5c65c9d..39a3b5e 100755 --- a/examples/structural-inference/tension/deterministic/optimize.py +++ b/examples/structural-inference/tension/deterministic/optimize.py @@ -43,7 +43,7 @@ # Don't try to optimize for the Young's modulus def make(n, eta, s0, R, d, **kwargs): """ - Maker with the Young's modulus fixed + Maker with the Young's modulus fixed """ return make_model(torch.tensor(0.5), n, eta, s0, R, d, device=device, **kwargs).to( device @@ -77,8 +77,10 @@ def make(n, eta, s0, R, d, **kwargs): # 3) Create the actual model model = optimize.DeterministicModel( - lambda *args, **kwargs: make(*args, block_size = time_chunk_size, **kwargs), - names, ics) + lambda *args, **kwargs: make(*args, block_size=time_chunk_size, **kwargs), + names, + ics, + ) # 4) Setup the optimizer niter = 5 diff --git a/examples/structural-inference/tension/statistical/demo_plot.py b/examples/structural-inference/tension/statistical/demo_plot.py index 6805d18..6605372 100755 --- a/examples/structural-inference/tension/statistical/demo_plot.py +++ b/examples/structural-inference/tension/statistical/demo_plot.py @@ -32,7 +32,7 @@ # Don't try to optimize for the Young's modulus def make(n, eta, s0, R, d, **kwargs): """ - Maker with the Young's modulus fixed + Maker with the Young's modulus fixed """ return make_model(torch.tensor(0.5), n, eta, s0, R, d, device=device, **kwargs).to( device @@ -57,7 +57,7 @@ def make(n, eta, s0, R, d, **kwargs): names = ["n", "eta", "s0", "R", "d"] sampler = optimize.StatisticalModel( - lambda *args, **kwargs: make(*args, block_size = time_chunk_size, **kwargs), + lambda *args, **kwargs: make(*args, block_size=time_chunk_size, **kwargs), names, [0.50, 0.49, 0.49, 0.48, 0.48], [0.02, 0.02, 0.03, 0.05, 0.05], diff --git a/examples/structural-inference/tension/statistical/infer.py b/examples/structural-inference/tension/statistical/infer.py index a998c41..1afa1eb 100755 --- a/examples/structural-inference/tension/statistical/infer.py +++ b/examples/structural-inference/tension/statistical/infer.py @@ -44,7 +44,7 @@ # Don't try to optimize for the Young's modulus def make(n, eta, s0, R, d, **kwargs): """ - Maker with Young's modulus fixed + Maker with Young's modulus fixed """ return make_model(torch.tensor(0.5), n, eta, s0, R, d, device=device, **kwargs).to( device @@ -89,9 +89,18 @@ def make(n, eta, s0, R, d, **kwargs): # 3) Create the actual model model = optimize.HierarchicalStatisticalModel( - lambda *args, **kwargs: make(*args, block_size = time_chunk_size, direct_solve_method = linear_solve_method, - **kwargs), - names, loc_loc_priors, loc_scale_priors, scale_scale_priors, eps, use_cached_guess = False + lambda *args, **kwargs: make( + *args, + block_size=time_chunk_size, + direct_solve_method=linear_solve_method, + **kwargs + ), + names, + loc_loc_priors, + loc_scale_priors, + scale_scale_priors, + eps, + use_cached_guess=False, ).to(device) # 4) Get the guide diff --git a/examples/structural-material-models/damage.py b/examples/structural-material-models/damage.py index 376d2e3..7da58a2 100755 --- a/examples/structural-material-models/damage.py +++ b/examples/structural-material-models/damage.py @@ -35,21 +35,20 @@ s0 = CP(torch.tensor(0.0)) theta0 = CP(torch.tensor(3000.0)) tau = CP(torch.tensor(100.0)) - + isotropic = hardening.Theta0VoceIsotropicHardeningModel(tau, theta0) kinematic = hardening.NoKinematicHardeningModel() - flowrule = flowrules.IsoKinViscoplasticity(n, eta, - s0, isotropic, kinematic) - + flowrule = flowrules.IsoKinViscoplasticity(n, eta, s0, isotropic, kinematic) + C = CP(torch.tensor(15.94)) A = CP(torch.tensor(-5379.0)) B = CP(torch.tensor(29316.3)) dmodel = damage.LarsonMillerDamage(C, A, B) - model = models.DamagedInelasticModel(E, flowrule, dmodel = dmodel) + model = models.DamagedInelasticModel(E, flowrule, dmodel=dmodel) - integrator = models.ModelIntegrator(model, block_size = 20) + integrator = models.ModelIntegrator(model, block_size=20) # Creep target_temperature = 550.0 + 273.15 @@ -74,7 +73,7 @@ state = integrator.solve_stress(times, stresses, temperatures) - strains = state[...,0] + strains = state[..., 0] plt.plot(times, strains) plt.xlabel("Time (hrs)") @@ -86,34 +85,40 @@ plt.ylabel("Stress (MPa)") plt.show() - plt.plot(times, state[...,-1]) + plt.plot(times, state[..., -1]) plt.xlabel("Time (hrs)") plt.ylabel("Damage") plt.show() - - # Simpler + + # Simpler R = CP(torch.tensor(1.0e-3)) dmodel = damage.ConstantDamage(R) - model = models.DamagedInelasticModel(E, flowrule, dmodel = dmodel) - integrator = models.ModelIntegrator(model, block_size = 20) + model = models.DamagedInelasticModel(E, flowrule, dmodel=dmodel) + integrator = models.ModelIntegrator(model, block_size=20) # Tension with unload - times = torch.cat([torch.linspace(0, 750, nsteps_hold) , torch.linspace(750, 751, nsteps_hold)[1:]]).unsqueeze(-1) - strains = torch.cat([torch.linspace(0, 0.4, nsteps_hold), torch.linspace(0.4, 0.398, nsteps_hold)[1:]]).unsqueeze(-1) + times = torch.cat( + [torch.linspace(0, 750, nsteps_hold), torch.linspace(750, 751, nsteps_hold)[1:]] + ).unsqueeze(-1) + strains = torch.cat( + [ + torch.linspace(0, 0.4, nsteps_hold), + torch.linspace(0.4, 0.398, nsteps_hold)[1:], + ] + ).unsqueeze(-1) temperatures = torch.zeros_like(strains) state = integrator.solve_strain(times, strains, temperatures) - stress = state[...,0] + stress = state[..., 0] - E_final = (stress[-2,0] - stress[-1,0])/(strains[-2,0] - strains[-1,0]) - d = state[-1,0,-1] + E_final = (stress[-2, 0] - stress[-1, 0]) / (strains[-2, 0] - strains[-1, 0]) + d = state[-1, 0, -1] print("Initial modulus: %f" % E.pvalue.cpu()) print("Final modulus: %f" % E_final.cpu()) - print("Calculated modulus: %f" % (E.pvalue * (1-d)).cpu()) + print("Calculated modulus: %f" % (E.pvalue * (1 - d)).cpu()) plt.plot(strains, stress) plt.xlabel("Strain (mm/mm)") plt.ylabel("Stress (MPa)") plt.show() - diff --git a/examples/structural-material-models/kocks-mecking-differentiable.py b/examples/structural-material-models/kocks-mecking-differentiable.py index 722a9a3..b0d15f4 100755 --- a/examples/structural-material-models/kocks-mecking-differentiable.py +++ b/examples/structural-material-models/kocks-mecking-differentiable.py @@ -74,23 +74,21 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): iso_hardening, kin_hardening, ) - - n_constant = CP(torch.tensor(15.0)) - eta_constant = CP(torch.tensor(1.0)) - ri_flowrule = flowrules.IsoKinViscoplasticity( - n_constant, eta_constant, s0, iso_hardening, kin_hardening + + ri_flowrule = flowrules.IsoKinRateIndependentPlasticity( + E, s0, iso_hardening, kin_hardening, s=10.0 ) - - sf = torch.tensor(10.0) + + sf = torch.tensor(50.0) flowrule = flowrules.SoftKocksMeckingRegimeFlowRule( ri_flowrule, rd_flowrule, g0, mu, b, eps0, k, sf ) model = models.InelasticModel(E, flowrule) - integrator = models.ModelIntegrator(model) + integrator = models.ModelIntegrator(model, block_size=50, linesearch=True) ngrid = 10 - nsteps = 200 + nsteps = 2000 elimits = torch.ones(ngrid) * 0.05 # Constant temperature, varying flow rate @@ -116,9 +114,18 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): ) results_temp = integrator.solve_strain(times, strains, temperatures) + + plt.figure() + plt.plot(strains, results_temp[:, :, 0], label=["%i K" % T for T in temps]) + plt.xlabel("Strain (mm/mm)") + plt.ylabel("Stress (MPa)") + plt.legend(loc="best") + plt.show() + yield_temp = calculate_yield(strains, results_temp[:, :, 0]) norm_temp = yield_temp / mu.value(temps) + plt.figure() plt.semilogy( g_rate.numpy(), norm_rate.numpy(), @@ -136,17 +143,11 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): label="Varying temperature", ) - grange = np.linspace(0.4,1.4,50) + grange = np.linspace(0.4, 1.4, 50) + plt.semilogy(grange, np.exp(C) * np.ones_like(grange), color="k", label=None) plt.semilogy( - grange, - np.exp(C)*np.ones_like(grange), - color = 'k', - label = None) - plt.semilogy( - grange, - np.exp(B.numpy()) * np.exp(A.numpy()*grange), - color = 'k', - label = None) + grange, np.exp(B.numpy()) * np.exp(A.numpy() * grange), color="k", label=None + ) plt.axvline(x=g0, ls="--", color="k") plt.legend(loc="best") @@ -154,3 +155,27 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): plt.ylabel("Normalized flow stress") plt.tight_layout() plt.show() + + plt.figure() + plt.plot(temperatures[0], yield_temp, "kx") + plt.xlabel("Temperature (K)") + plt.ylabel("Yield stress (MPa)") + plt.show() + + nstress = 5 + times, stresses, temperatures, cycles = experiments.make_creep_tests( + torch.ones_like(temps) * 150.0, + temps, + torch.ones_like(temps), + torch.ones_like(temps) * 1000.0, + 50, + 150, + ) + + creep_strains = integrator.solve_stress(times, stresses, temperatures)[:, :, 0] + plt.figure() + plt.semilogy(times, creep_strains, label=["%i K" % T for T in temps]) + plt.xlabel("Time (s)") + plt.ylabel("Creep strain (mm/mm)") + plt.legend(loc="best") + plt.show() diff --git a/examples/structural-material-models/kocks-mecking.py b/examples/structural-material-models/kocks-mecking.py index 3b61444..26aec92 100755 --- a/examples/structural-material-models/kocks-mecking.py +++ b/examples/structural-material-models/kocks-mecking.py @@ -74,7 +74,7 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): iso_hardening, kin_hardening, ) - + n_constant = CP(torch.tensor(15.0)) eta_constant = CP(torch.tensor(1.0)) ri_flowrule = flowrules.IsoKinViscoplasticity( @@ -135,17 +135,11 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): label="Varying temperature", ) - grange = np.linspace(0.4,1.4,50) - plt.semilogy( - grange, - np.exp(C)*np.ones_like(grange), - color = 'k', - label = None) + grange = np.linspace(0.4, 1.4, 50) + plt.semilogy(grange, np.exp(C) * np.ones_like(grange), color="k", label=None) plt.semilogy( - grange, - np.exp(B.numpy()) * np.exp(A.numpy()*grange), - color = 'k', - label = None) + grange, np.exp(B.numpy()) * np.exp(A.numpy() * grange), color="k", label=None + ) plt.axvline(x=g0, ls="--", color="k") plt.legend(loc="best") diff --git a/examples/structural-material-models/rate_independent.py b/examples/structural-material-models/rate_independent.py index 6057a4c..2a4c18c 100755 --- a/examples/structural-material-models/rate_independent.py +++ b/examples/structural-material-models/rate_independent.py @@ -13,6 +13,7 @@ import torch import numpy as np +from matplotlib.lines import Line2D import matplotlib.pyplot as plt from pyoptmat import models, flowrules, experiments, hardening @@ -23,7 +24,7 @@ if __name__ == "__main__": nrates = 5 - erates = torch.logspace(-1, -5, nrates) + erates = torch.logspace(-1, -8, nrates) temps = torch.ones_like(erates) elimits = torch.ones_like(erates) * 0.1 @@ -36,7 +37,7 @@ # Viscoplastic base model E = CP(torch.tensor(100000.0)) - n = CP(torch.tensor(5.0)) + n = CP(torch.tensor(15.0)) eta = CP(torch.tensor(50.0)) s0 = CP(torch.tensor(25.0)) @@ -50,19 +51,17 @@ n, eta, s0, iso_hardening, kin_hardening ) - # Approximately rate independent flow rule - lmbda = torch.tensor(0.999) - eps_ref = torch.tensor(1e-10) - - flowrule = flowrules.RateIndependentFlowRuleWrapper(base_flowrule, lmbda, eps_ref) + flowrule = flowrules.IsoKinRateIndependentPlasticity( + E, s0, iso_hardening, kin_hardening, s=10.0 + ) base_model = models.InelasticModel(E, base_flowrule) - base_integrator = models.ModelIntegrator(base_model) + base_integrator = models.ModelIntegrator(base_model, block_size=10, linesearch=True) rd_stresses = base_integrator.solve_strain(times, strains, temperatures)[:, :, 0] model = models.InelasticModel(E, flowrule) - integrator = models.ModelIntegrator(model) + integrator = models.ModelIntegrator(model, block_size=10, linesearch=True) ri_stresses = integrator.solve_strain(times, strains, temperatures)[:, :, 0] @@ -81,12 +80,14 @@ def analytic(strain): return pred for i in range(nrates): - plt.plot(strains[:, i], ri_stresses[:, i], "k-", lw=3) - plt.plot(strains[:, i], rd_stresses[:, i], "k--", lw=3) + l1 = plt.plot(strains[:, i], ri_stresses[:, i], "k-", lw=3) + l2 = plt.plot(strains[:, i], rd_stresses[:, i], "k--", lw=3) enp = strains[:, i].numpy() exact = analytic(enp) - plt.plot(enp, exact, color="r", lw=2, ls=":") + l3 = plt.plot(enp, exact, color="r", lw=2, ls=":") + + plt.gca().legend([l1[0], l2[0], l3[0]], ["Approx RI", "RD", "Analytic"], loc="best") plt.xlabel("Strain (mm/mm)") plt.ylabel("Stress (MPa)") diff --git a/examples/structural-material-models/rate_independent_approximations.py b/examples/structural-material-models/rate_independent_approximations.py deleted file mode 100755 index 7753dd5..0000000 --- a/examples/structural-material-models/rate_independent_approximations.py +++ /dev/null @@ -1,125 +0,0 @@ -#!/usr/bin/env python3 - -""" - This example demonstrates the ability of pyoptmat to integrate material models under - strain control, stress control, or a combination of both. - - The example first simulates standard, strain-rate controlled tensile curves - for a simple material model for Alloy 617 developed by Messner, Phan, and Sham. These - simulations then take time, strain, and temperature as inputs and outputs the - resulting stresses. - - Then, the example uses the stresses from the strain controlled simulations - to repeat the simulations under stress control. These simulations therefore - take as input time, temperature, and stress and output strain. The plot demonstrates that - the strain controlled and stress controlled simulations produce identical results. - - Finally, the example integrates the model twice for each temperature, once in strain - control and once in stress control, but integrated all the experiments at once. - This example then demos the ability of pyoptmat to integrate both stress and strain - controlled tests in the same vectorized calculation. -""" - -import sys - -sys.path.append("../..") - -import torch - -import matplotlib.pyplot as plt - -from pyoptmat import models, flowrules, temperature, experiments, hardening - -torch.set_default_tensor_type(torch.DoubleTensor) - -if __name__ == "__main__": - # This example illustrates the temperature and rate sensitivity of - # Alloy 617 with the model given in: - # Messner, Phan, and Sham. IJPVP 2019. - - # Setup test conditions - ntemps = 5 - temps = torch.linspace(750, 950, ntemps) + 273.15 - elimits = torch.ones(ntemps) * 0.25 - erates = torch.logspace(-16,-1,ntemps) - - nsteps = 200 - - times, strains, temperatures, cycles = experiments.make_tension_tests( - erates, temps, elimits, nsteps - ) - - # Perfect viscoplastic model - C = torch.tensor(-5.0) - mu = temperature.PolynomialScaling( - torch.tensor([-1.34689305e-02, -5.18806776e00, 7.86708330e04]) - ) - E = temperature.PolynomialScaling( - torch.tensor([-3.48056033e-02, -1.44398964e01, 2.06464967e05]) - ) - sy = temperature.ConstantParameter(200.0) - - C = torch.tensor(12000.0) - g = torch.tensor(10.1) - kh = hardening.FAKinematicHardeningModel(temperature.ConstantParameter(C), temperature.ConstantParameter(g)) - - ih = hardening.Theta0VoceIsotropicHardeningModel(temperature.ConstantParameter(100.0), temperature.ConstantParameter(500.0)) - - fr = flowrules.IsoKinRateIndependentPlasticity( - E, sy, ih, kh - ) - model = models.InelasticModel(E, fr) - integrator = models.ModelIntegrator(model, guess_type = "previous", linesearch = True, block_size = 10) - - stresses = integrator.solve_strain(times, strains, temperatures)[:, :, 0] - - # Now do it backwards! - strains_prime = integrator.solve_stress(times, stresses, temperatures)[:, :, 0] - - plt.figure() - for ei, epi, Ti, si in zip( - strains.T.numpy(), - strains_prime.T.numpy(), - temperatures.T.numpy(), - stresses.T.numpy(), - ): - (l,) = plt.plot(ei, si, label="T = %3.0fK" % Ti[0]) - plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw = 4, alpha = 0.5) - - plt.legend(loc="best") - plt.xlabel("Strain (mm/mm)") - plt.ylabel("Stress (Mpa)") - plt.title("Comparison between strain and stress control") - plt.show() - - # Now do it both ways at once! - big_times = torch.cat((times, times), 1) - big_temperatures = torch.cat((temperatures, temperatures), 1) - big_data = torch.cat((strains, stresses), 1) - big_types = torch.zeros(2 * ntemps) - big_types[ntemps:] = 1 - - both_results = integrator.solve_both( - big_times, big_temperatures, big_data, big_types - ) - - plt.figure() - for i in range(ntemps): - (l,) = plt.plot( - strains[:, i], both_results[:, i, 0], label="T = %3.0fK" % Ti[i] - ) - plt.plot( - both_results[:, i + ntemps, 0], - stresses[:, i], - label=None, - ls="--", - color=l.get_color(), - lw = 4, - alpha = 0.5 - ) - - plt.xlabel("Strain (mm/mm)") - plt.ylabel("Stress (MPa)") - plt.legend(loc="best") - plt.title("Running both strain and stress control at once") - plt.show() diff --git a/examples/structural-material-models/rate_independent_stability.py b/examples/structural-material-models/rate_independent_stability.py index 4abca52..cb6cd7e 100755 --- a/examples/structural-material-models/rate_independent_stability.py +++ b/examples/structural-material-models/rate_independent_stability.py @@ -21,11 +21,13 @@ torch.set_default_tensor_type(torch.DoubleTensor) if __name__ == "__main__": - nrates = 5 + # Simple stability test for rate independent model - erates = torch.logspace(-1, -6, nrates) - temps = torch.ones_like(erates) - elimits = torch.ones_like(erates) * 0.1 + # Setup test conditions + nrates = 5 + elimits = torch.ones(nrates) * 0.25 + erates = torch.logspace(-8, -1, nrates) + temps = np.ones_like(erates) nsteps = 200 @@ -36,58 +38,91 @@ # Viscoplastic base model E = CP(torch.tensor(100000.0)) - n = CP(torch.tensor(5.0)) - eta = CP(torch.tensor(50.0)) - s0 = CP(torch.tensor(25.0)) - + s0 = CP(torch.tensor(100.0)) R = CP(torch.tensor(100.0)) d = CP(torch.tensor(20.0)) iso_hardening = hardening.VoceIsotropicHardeningModel(R, d) kin_hardening = hardening.NoKinematicHardeningModel() - - base_flowrule = flowrules.IsoKinViscoplasticity( - n, eta, s0, iso_hardening, kin_hardening + fr = flowrules.IsoKinRateIndependentPlasticity( + E, s0, iso_hardening, kin_hardening, s=0.5 ) + model = models.InelasticModel(E, fr) + integrator = models.ModelIntegrator(model, linesearch=True, block_size=20) + + stresses = integrator.solve_strain(times, strains, temperatures)[:, :, 0] + + # Now do it backwards! + strains_prime = integrator.solve_stress(times, stresses, temperatures)[:, :, 0] + + plt.figure() + for i, (ei, epi, Ti, si) in enumerate( + zip( + strains.T.numpy(), + strains_prime.T.numpy(), + temperatures.T.numpy(), + stresses.T.numpy(), + ) + ): + (l,) = plt.plot(ei, si, label=r"$\dot{\varepsilon}$ = %3.0e 1/s" % erates[i]) + plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw=4, alpha=0.5) + + plt.legend(loc="best") + plt.xlabel("Strain (mm/mm)") + plt.ylabel("Stress (Mpa)") + plt.title("Comparison between strain and stress control") + plt.show() - # Approximately rate independent flow rule - n_ri = CP(torch.tensor(15)) - eta_ri = CP(torch.tensor(1.0)) - flowrule = flowrules.IsoKinViscoplasticity(n_ri, eta_ri, s0, - iso_hardening, kin_hardening) - - base_model = models.InelasticModel(E, base_flowrule) - base_integrator = models.ModelIntegrator(base_model) - - model = models.InelasticModel(E, flowrule) - integrator = models.ModelIntegrator(model) - - ri_stresses = integrator.solve_strain(times, strains, temperatures)[:, :, 0] - ri_strains = integrator.solve_stress(times, ri_stresses, temperatures)[:, :, 0] - - def analytic(strain): - Ev = E.value(0).numpy() - s0v = s0.value(0).numpy() - Rv = R.value(0).numpy() - dv = d.value(0).numpy() - - pred = Ev * strain - plastic = pred > s0v - ystrain = s0v / Ev - pstrain = strain[plastic] - ystrain - pred[plastic] = s0v + Rv * (1.0 - np.exp(-dv * pstrain)) + # Now do it both ways at once! + big_times = torch.cat((times, times), 1) + big_temperatures = torch.cat((temperatures, temperatures), 1) + big_data = torch.cat((strains, stresses), 1) + big_types = torch.zeros(2 * nrates) + big_types[nrates:] = 1 - return pred + both_results = integrator.solve_both( + big_times, big_temperatures, big_data, big_types + ) + plt.figure() for i in range(nrates): - plt.plot(strains[:, i], ri_stresses[:, i], "k-", lw=3) - plt.plot(ri_strains[:, i], ri_stresses[:, i], "b-", lw = 2) - - enp = strains[:, i].numpy() - exact = analytic(enp) - plt.plot(enp, exact, color="r", lw=2, ls=":") + (l,) = plt.plot( + strains[:, i], + both_results[:, i, 0], + label=r"$\dot{\varepsilon}$ = %3.0e 1/s" % erates[i], + ) + plt.plot( + both_results[:, i + nrates, 0], + stresses[:, i], + label=None, + ls="--", + color=l.get_color(), + lw=4, + alpha=0.5, + ) plt.xlabel("Strain (mm/mm)") plt.ylabel("Stress (MPa)") + plt.legend(loc="best") + plt.title("Running both strain and stress control at once") + plt.show() + + nstress = 5 + input_stress = torch.linspace(50, 150, nstress) + times, stresses, temperatures, cycles = experiments.make_creep_tests( + input_stress, + torch.ones_like(input_stress), + torch.ones_like(input_stress), + torch.ones_like(input_stress) * 1000.0, + 50, + 50, + ) + + creep_strains = integrator.solve_stress(times, stresses, temperatures)[:, :, 0] + plt.figure() + plt.plot(times, creep_strains, label=["%i MPa" % s for s in input_stress]) + plt.xlabel("Time (s)") + plt.ylabel("Creep strain (mm/mm)") + plt.legend(loc="best") plt.show() diff --git a/examples/structural-material-models/stress_strain_control.py b/examples/structural-material-models/stress_strain_control.py index a4234e3..39c66e6 100755 --- a/examples/structural-material-models/stress_strain_control.py +++ b/examples/structural-material-models/stress_strain_control.py @@ -91,7 +91,7 @@ stresses.T.numpy(), ): (l,) = plt.plot(ei, si, label="T = %3.0fK" % Ti[0]) - plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw = 4, alpha = 0.5) + plt.plot(epi, si, label=None, ls="--", color=l.get_color(), lw=4, alpha=0.5) plt.legend(loc="best") plt.xlabel("Strain (mm/mm)") @@ -120,7 +120,9 @@ stresses[:, i], label=None, ls="--", - color=l.get_color(), lw = 4, alpha = 0.5 + color=l.get_color(), + lw=4, + alpha=0.5, ) plt.xlabel("Strain (mm/mm)") diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index f9f3ec4..43150d1 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -20,8 +20,14 @@ def newton_raphson_chunk( - fn, x0, solver, rtol=1e-6, atol=1e-10, miter=100, throw_on_fail=False, - linesearch = False + fn, + x0, + solver, + rtol=1e-6, + atol=1e-10, + miter=100, + throw_on_fail=False, + linesearch=False, ): """ Solve a nonlinear system with Newton's method with a tensor for a @@ -49,7 +55,9 @@ def newton_raphson_chunk( nR0 = nR i = 0 - while (i < miter) and torch.any(torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol))): + while (i < miter) and torch.any( + torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol)) + ): print(i, torch.max(nR)) dx = solver.solve(J, R) if linesearch: @@ -57,9 +65,9 @@ def newton_raphson_chunk( else: x -= dx R, J = fn(x) - nR = torch.norm(R, dim = -1) + nR = torch.norm(R, dim=-1) i += 1 - + if i == miter: if throw_on_fail: raise RuntimeError("Implicit solve did not succeed.") @@ -67,11 +75,14 @@ def newton_raphson_chunk( return x -def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=0.0, miter = 10): + +def chunk_linesearch( + x, dx, fn, R0, overall_rtol, overall_atol, sigma=2.0, c=0.0, miter=10 +): """ Backtracking linesearch for the chunk NR algorithm. - Terminates when the Armijo criteria is reached, or you exceed some maximum iterations, or when you would meet the + Terminates when the Armijo criteria is reached, or you exceed some maximum iterations, or when you would meet the convergence requirements with the current alpha Args: @@ -86,20 +97,24 @@ def chunk_linesearch(x, dx, fn, R0, overall_rtol, overall_atol, sigma = 2.0, c=0 c (scalar): stopping criteria miter (scalar): maximum iterations """ - alpha = torch.ones(x.shape[:-1], device = x.device) - nR0 = torch.norm(R0, dim = -1) + alpha = torch.ones(x.shape[:-1], device=x.device) + nR0 = torch.norm(R0, dim=-1) i = 0 while True: R, J = fn(x - dx * alpha.unsqueeze(-1)) - nR = torch.norm(R, dim = -1)**2.0 - crit = nR0**2.0 + 2.0 * c * alpha * torch.einsum('...i,...i', R0, dx) + nR = torch.norm(R, dim=-1) ** 2.0 + crit = nR0**2.0 + 2.0 * c * alpha * torch.einsum("...i,...i", R0, dx) i += 1 - if torch.all(nR < crit) or i >= miter or torch.all(torch.sqrt(nR) < overall_atol) or torch.all(torch.sqrt(nR)/nR0 < overall_rtol): + if ( + torch.all(nR < crit) + or i >= miter + or torch.all(torch.sqrt(nR) < overall_atol) + or torch.all(torch.sqrt(nR) / nR0 < overall_rtol) + ): break alpha[nR >= crit] /= sigma - return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim = -1), alpha - + return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim=-1), alpha class BidiagonalOperator(torch.nn.Module): diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index 595d615..6073623 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -1125,9 +1125,10 @@ def history_rate(self, s, h, t, T, e): """ return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) + class PerfectRateIndependentPlasticity(FlowRule): """ - Perfect rate independent plasticity + Perfect rate independent plasticity Args: sy (|TP|): yield stress @@ -1171,7 +1172,7 @@ def flow_rate(self, s, h, t, T, e): fr = torch.zeros_like(s) yielding = self.f(s, h, T) > 0 fr[yielding] = e[yielding] - + dfr = torch.zeros_like(s) return fr, dfr @@ -1227,6 +1228,7 @@ def history_rate(self, s, h, t, T, e): """ return torch.zeros_like(h), torch.zeros(h.shape + h.shape[-1:]) + class IsoKinViscoplasticity(FlowRule): """ Viscoplasticity with isotropic and kinematic hardening, defined as @@ -1530,8 +1532,8 @@ def dhist_derate(self, s, h, t, T, e): class IsoKinRateIndependentPlasticity(FlowRule): """ - An approximation to rate independent plasticity. The model is defined by - + An approximation to rate independent plasticity. The model is defined by + .. math:: \\dot{\\varepsilon}_{in} = \\xi(f) \\dot{\\varepsilon}_{p,ri} @@ -1571,9 +1573,9 @@ class IsoKinRateIndependentPlasticity(FlowRule): s (float): scale factor for the sigmoid function, controls the amount of smoothing at the onset of plasticity """ - def __init__(self, E, sy, isotropic, kinematic, soffset = 1e-10, s = 1.0): + def __init__(self, E, sy, isotropic, kinematic, soffset=1e-10, s=1.0): super().__init__() - self. E = E + self.E = E self.isotropic = isotropic self.kinematic = kinematic self.sy = sy @@ -1595,7 +1597,7 @@ def f(self, s, h, T): """ ih = self.isotropic.value(h[..., : self.isotropic.nhist]) kh = self.kinematic.value(h[..., self.isotropic.nhist :]) - return torch.abs(s-kh) - self.sy(T) - ih + return torch.abs(s - kh) - self.sy(T) - ih def df_ds(self, s, h, T): """ @@ -1610,7 +1612,7 @@ def df_ds(self, s, h, T): the derivative value """ kh = self.kinematic.value(h[..., self.isotropic.nhist :]) - return torch.sign(s-kh) + return torch.sign(s - kh) def df_dh(self, s, h, T): """ @@ -1626,9 +1628,11 @@ def df_dh(self, s, h, T): """ kh = self.kinematic.value(h[..., self.isotropic.nhist :]) di = -self.isotropic.dvalue(h[..., : self.isotropic.nhist]) - dk = -torch.sign(s-kh).unsqueeze(-1) * self.kinematic.dvalue(h[..., self.isotropic.nhist :]) + dk = -torch.sign(s - kh).unsqueeze(-1) * self.kinematic.dvalue( + h[..., self.isotropic.nhist :] + ) - return torch.cat([di,dk], axis = -1) + return torch.cat([di, dk], axis=-1) def ep_residual(self, ep, s, h, t, T, e): """ @@ -1651,10 +1655,28 @@ def ep_residual(self, ep, s, h, t, T, e): kh = self.kinematic.value(hkin) - i_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.history_rate(s, hiso, t, ep, T, e).unsqueeze(-1)).squeeze(-1).squeeze(-1) - k_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.history_rate(s, hkin, t, ep, T, e).unsqueeze(-1)).squeeze(-1).squeeze(-1) + i_dot = ( + utility.mbmm( + self.isotropic.dvalue(hiso).unsqueeze(-2), + self.isotropic.history_rate(s, hiso, t, ep, T, e).unsqueeze(-1), + ) + .squeeze(-1) + .squeeze(-1) + ) + k_dot = ( + utility.mbmm( + self.kinematic.dvalue(hkin).unsqueeze(-2), + self.kinematic.history_rate(s, hkin, t, ep, T, e).unsqueeze(-1), + ) + .squeeze(-1) + .squeeze(-1) + ) - return torch.sign(s - kh + self.soffset) * self.E(T) * (e - ep) - torch.sign(s - kh + self.soffset) * k_dot - i_dot + return ( + torch.sign(s - kh + self.soffset) * self.E(T) * (e - ep) + - torch.sign(s - kh + self.soffset) * k_dot + - i_dot + ) def ep_jacobian_ep(self, ep, s, h, t, T, e): """ @@ -1675,11 +1697,21 @@ def ep_jacobian_ep(self, ep, s, h, t, T, e): hkin = h[..., self.isotropic.nhist :] kh = self.kinematic.value(hkin) - - di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_derate(s, hiso, t, ep, T, e)).squeeze(-1) - dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_derate(s, hkin, t, ep, T, e)).squeeze(-1) - return -(torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) - torch.sign(s - kh + self.soffset).unsqueeze(-1) * dk_dot - di_dot + di_dot = utility.mbmm( + self.isotropic.dvalue(hiso).unsqueeze(-2), + self.isotropic.dhistory_rate_derate(s, hiso, t, ep, T, e), + ).squeeze(-1) + dk_dot = utility.mbmm( + self.kinematic.dvalue(hkin).unsqueeze(-2), + self.kinematic.dhistory_rate_derate(s, hkin, t, ep, T, e), + ).squeeze(-1) + + return ( + -(torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) + - torch.sign(s - kh + self.soffset).unsqueeze(-1) * dk_dot + - di_dot + ) def ep_jacobian_s(self, ep, s, h, t, T, e): """ @@ -1700,9 +1732,15 @@ def ep_jacobian_s(self, ep, s, h, t, T, e): hkin = h[..., self.isotropic.nhist :] kh = self.kinematic.value(hkin) - - di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_dstress(s, hiso, t, ep, T, e).unsqueeze(-1)).squeeze(-1) - dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_dstress(s, hkin, t, ep, T, e).unsqueeze(-1)).squeeze(-1) + + di_dot = utility.mbmm( + self.isotropic.dvalue(hiso).unsqueeze(-2), + self.isotropic.dhistory_rate_dstress(s, hiso, t, ep, T, e).unsqueeze(-1), + ).squeeze(-1) + dk_dot = utility.mbmm( + self.kinematic.dvalue(hkin).unsqueeze(-2), + self.kinematic.dhistory_rate_dstress(s, hkin, t, ep, T, e).unsqueeze(-1), + ).squeeze(-1) return -torch.sign(s - kh + self.soffset).unsqueeze(-1) * dk_dot - di_dot @@ -1726,10 +1764,16 @@ def ep_jacobian_h(self, ep, s, h, t, T, e): kh = self.kinematic.value(hkin) - di_dot = utility.mbmm(self.isotropic.dvalue(hiso).unsqueeze(-2), self.isotropic.dhistory_rate_dhistory(s, hiso, t, ep, T, e)).squeeze(-2) - dk_dot = utility.mbmm(self.kinematic.dvalue(hkin).unsqueeze(-2), self.kinematic.dhistory_rate_dhistory(s, hkin, t, ep, T, e)).squeeze(-2) + di_dot = utility.mbmm( + self.isotropic.dvalue(hiso).unsqueeze(-2), + self.isotropic.dhistory_rate_dhistory(s, hiso, t, ep, T, e), + ).squeeze(-2) + dk_dot = utility.mbmm( + self.kinematic.dvalue(hkin).unsqueeze(-2), + self.kinematic.dhistory_rate_dhistory(s, hkin, t, ep, T, e), + ).squeeze(-2) - return torch.cat([-di_dot, -torch.sign(s - kh).unsqueeze(-1) * dk_dot], dim = -1) + return torch.cat([-di_dot, -torch.sign(s - kh).unsqueeze(-1) * dk_dot], dim=-1) def ep_jacobian_e(self, ep, s, h, t, T, e): """ @@ -1749,7 +1793,7 @@ def ep_jacobian_e(self, ep, s, h, t, T, e): hkin = h[..., self.isotropic.nhist :] kh = self.kinematic.value(hkin) - + return (torch.sign(s - kh + self.soffset) * self.E(T)).unsqueeze(-1) def sig(self, x): @@ -1765,7 +1809,7 @@ def sig(self, x): Args: x (torch.tensor): input to sigmoid """ - return (torch.tanh(self.s * x) + 1.0)/2.0 + return (torch.tanh(self.s * x) + 1.0) / 2.0 def dsig(self, x): """ @@ -1774,7 +1818,7 @@ def dsig(self, x): Args: x (torch.tensor): input to sigmoid """ - return 0.5*self.s / torch.cosh(self.s * x)**2.0 + return 0.5 * self.s / torch.cosh(self.s * x) ** 2.0 def mix_fn(self, s, h, T): """ @@ -1820,8 +1864,11 @@ def plastic_rate(self, s, h, t, T, e): T (torch.tensor): temperature e (torch.tensor): total strain rate """ + def RJ(x): - return self.ep_residual(x, s+self.soffset, h, t, T, e), self.ep_jacobian_ep(x, s+self.soffset, h, t, T, e).squeeze(-1) + return self.ep_residual( + x, s + self.soffset, h, t, T, e + ), self.ep_jacobian_ep(x, s + self.soffset, h, t, T, e).squeeze(-1) return solvers.scalar_newton(RJ, torch.zeros_like(e)) @@ -1845,10 +1892,16 @@ def flow_rate(self, s, h, t, T, e): bfr = self.plastic_rate(s, h, t, T, e) fr = mf * bfr - dfr = -(self.ep_jacobian_s(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)).squeeze(-1) * mf + self.dmix_fn_ds(s, h, T) * bfr - + dfr = ( + -( + self.ep_jacobian_s(bfr, s, h, t, T, e) + / self.ep_jacobian_ep(bfr, s, h, t, T, e) + ).squeeze(-1) + * mf + + self.dmix_fn_ds(s, h, T) * bfr + ) + return fr, dfr - def dflow_derate(self, s, h, t, T, e): """ @@ -1870,7 +1923,13 @@ def dflow_derate(self, s, h, t, T, e): """ pfr = self.plastic_rate(s, h, t, T, e) mf = self.mix_fn(s, h, T) - dfr = -(self.ep_jacobian_e(pfr, s, h, t, T, e) / self.ep_jacobian_ep(pfr, s, h, t, T, e)).squeeze(-1) * mf + dfr = ( + -( + self.ep_jacobian_e(pfr, s, h, t, T, e) + / self.ep_jacobian_ep(pfr, s, h, t, T, e) + ).squeeze(-1) + * mf + ) return dfr @@ -1890,8 +1949,11 @@ def dflow_dhist(self, s, h, t, T, e): """ bfr = self.plastic_rate(s, h, t, T, e) mf = self.mix_fn(s, h, T) - dfr = -(self.ep_jacobian_h(bfr, s, h, t, T, e) / self.ep_jacobian_ep(bfr, s, h, t, T, e)) * mf.unsqueeze(-1) + self.dmix_fn_dh(s, h, T) * bfr.unsqueeze(-1) - + dfr = -( + self.ep_jacobian_h(bfr, s, h, t, T, e) + / self.ep_jacobian_ep(bfr, s, h, t, T, e) + ) * mf.unsqueeze(-1) + self.dmix_fn_dh(s, h, T) * bfr.unsqueeze(-1) + return dfr.unsqueeze(-2) @property @@ -1937,8 +1999,8 @@ def history_rate(self, s, h, t, T, e): ) df_dh = self.dflow_dhist(s, h, t, T, e) - df_di = df_dh[...,:self.isotropic.nhist] - df_dk = df_dh[...,self.isotropic.nhist:] + df_di = df_dh[..., : self.isotropic.nhist] + df_dk = df_dh[..., self.isotropic.nhist :] # Jacobian contribution hdiv = torch.cat( @@ -1950,7 +2012,7 @@ def history_rate(self, s, h, t, T, e): self.isotropic.dhistory_rate_derate( s, hiso, t, erate, T, e ), - df_di + df_di, ), self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e) * df_dk, @@ -1963,7 +2025,7 @@ def history_rate(self, s, h, t, T, e): self.kinematic.dhistory_rate_derate( s, hkin, t, erate, T, e ), - df_di + df_di, ), self.kinematic.dhistory_rate_dhistory(s, hkin, t, erate, T, e) + utility.mbmm( @@ -2038,15 +2100,16 @@ def dhist_derate(self, s, h, t, T, e): return torch.cat( [ - self.isotropic.dhistory_rate_dtotalrate(s, hiso, t, erate, T, e) + utility.mbmm( + self.isotropic.dhistory_rate_dtotalrate(s, hiso, t, erate, T, e) + + utility.mbmm( self.isotropic.dhistory_rate_derate(s, hiso, t, erate, T, e), derate[..., None, None], - )[...,0], - self.kinematic.dhistory_rate_dtotalrate(s, hkin, t, erate, T, e) + utility.mbmm( + )[..., 0], + self.kinematic.dhistory_rate_dtotalrate(s, hkin, t, erate, T, e) + + utility.mbmm( self.kinematic.dhistory_rate_derate(s, hkin, t, erate, T, e), - derate[..., None, None])[...,0], + derate[..., None, None], + )[..., 0], ], dim=-1, ) - - diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 282e17e..13a0f93 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -133,7 +133,7 @@ def forward(self, t, y, erate, T): # Logically we should return the derivative wrt T, but right now # we're not going to use it Trate = torch.zeros_like(y) - + return result, dresult, drate, Trate @@ -283,7 +283,7 @@ class ModelIntegrator(nn.Module): """ - def __init__(self, model, *args, use_adjoint=True, bisect_first = False, **kwargs): + def __init__(self, model, *args, use_adjoint=True, bisect_first=False, **kwargs): super().__init__(*args) self.model = model self.use_adjoint = use_adjoint @@ -324,7 +324,7 @@ def solve_both(self, times, temperatures, idata, control): idata, temperatures, control, - bisect_first = self.bisect_first + bisect_first=self.bisect_first, ) return self.imethod(bmodel, init, times, **self.kwargs_for_integration) @@ -433,7 +433,18 @@ class BothBasedModel(nn.Module): indices: split into strain and stress control """ - def __init__(self, model, times, rates, base, temps, control, bisect_first = False, *args, **kwargs): + def __init__( + self, + model, + times, + rates, + base, + temps, + control, + bisect_first=False, + *args, + **kwargs + ): super().__init__(*args, **kwargs) self.model = model self.control = control @@ -441,15 +452,27 @@ def __init__(self, model, times, rates, base, temps, control, bisect_first = Fal self.econtrol = self.control == 0 self.scontrol = self.control == 1 - self.emodel = StrainBasedModel(self.model, - utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.econtrol], rates[...,self.econtrol]), - utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.econtrol], temps[...,self.econtrol])) + self.emodel = StrainBasedModel( + self.model, + utility.ArbitraryBatchTimeSeriesInterpolator( + times[..., self.econtrol], rates[..., self.econtrol] + ), + utility.ArbitraryBatchTimeSeriesInterpolator( + times[..., self.econtrol], temps[..., self.econtrol] + ), + ) self.smodel = StressBasedModel( - self.model, - utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], rates[...,self.scontrol]), - utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], base[...,self.scontrol]), - utility.ArbitraryBatchTimeSeriesInterpolator(times[...,self.scontrol], temps[...,self.scontrol]), - bisect_first = bisect_first + self.model, + utility.ArbitraryBatchTimeSeriesInterpolator( + times[..., self.scontrol], rates[..., self.scontrol] + ), + utility.ArbitraryBatchTimeSeriesInterpolator( + times[..., self.scontrol], base[..., self.scontrol] + ), + utility.ArbitraryBatchTimeSeriesInterpolator( + times[..., self.scontrol], temps[..., self.scontrol] + ), + bisect_first=bisect_first, ) def forward(self, t, y): @@ -464,18 +487,22 @@ def forward(self, t, y): n = (y.shape[-1],) base = y.shape[:-1] - actual_rates = torch.zeros(base + n, device = t.device) - actual_jacs = torch.zeros(base + n + n, device = t.device) + actual_rates = torch.zeros(base + n, device=t.device) + actual_jacs = torch.zeros(base + n + n, device=t.device) if torch.any(self.econtrol): - strain_rates, strain_jacs = self.emodel(t[...,self.econtrol], y[...,self.econtrol,:]) - actual_rates[...,self.econtrol,:] = strain_rates - actual_jacs[...,self.econtrol,:,:] = strain_jacs + strain_rates, strain_jacs = self.emodel( + t[..., self.econtrol], y[..., self.econtrol, :] + ) + actual_rates[..., self.econtrol, :] = strain_rates + actual_jacs[..., self.econtrol, :, :] = strain_jacs if torch.any(self.scontrol): - stress_rates, stress_jacs = self.smodel(t[...,self.scontrol], y[...,self.scontrol,:]) - actual_rates[...,self.scontrol,:] = stress_rates - actual_jacs[...,self.scontrol,:,:] = stress_jacs + stress_rates, stress_jacs = self.smodel( + t[..., self.scontrol], y[..., self.scontrol, :] + ) + actual_rates[..., self.scontrol, :] = stress_rates + actual_jacs[..., self.scontrol, :, :] = stress_jacs return actual_rates, actual_jacs @@ -519,7 +546,19 @@ class StressBasedModel(nn.Module): T_fn: T(t) """ - def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate = 1e3, bisect_first = False, *args, **kwargs): + def __init__( + self, + model, + srate_fn, + stress_fn, + T_fn, + min_erate=-1e2, + max_erate=1e3, + guess_erate=1.0e-3, + bisect_first=False, + *args, + **kwargs + ): super().__init__(*args, **kwargs) self.model = model self.srate_fn = srate_fn @@ -528,6 +567,7 @@ def __init__(self, model, srate_fn, stress_fn, T_fn, min_erate = -1e2, max_erate self.min_erate = min_erate self.max_erate = max_erate self.bisect_first = bisect_first + self.guess_erate = guess_erate def forward(self, t, y): """ @@ -542,37 +582,46 @@ def forward(self, t, y): cT = self.T_fn(t) def RJ(erate): - ydot, _, Je, _ = self.model(t, - torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), - erate, cT) + ydot, _, Je, _ = self.model( + t, torch.cat([cs.unsqueeze(-1), y[..., 1:]], dim=-1), erate, cT + ) R = ydot[..., 0] - csr J = Je[..., 0] return R, J - - # Doing the detach here actually makes the parameter gradients wrong... + + # Doing the detach here actually makes the parameter gradients wrong... if self.bisect_first: - erate = solvers.scalar_bisection_newton(RJ, - torch.ones_like(y[...,0]) * self.min_erate, - torch.ones_like(y[...,0]) * self.max_erate) + erate = solvers.scalar_bisection_newton( + RJ, + torch.ones_like(y[..., 0]) * self.min_erate, + torch.ones_like(y[..., 0]) * self.max_erate, + ) else: - erate = solvers.scalar_newton(RJ, - torch.zeros_like(y[...,0])) - - ydot, J, Je, _ = self.model(t, - torch.cat([cs.unsqueeze(-1), y[...,1:]], dim = -1), - erate, cT) - - # There is an annoying extra term that is the derivative of the history rate with respect to the - # solved for strain rate times the derivative of the strain rate with respect to history - t1 = Je[...,1:].unsqueeze(-1) - t2 = utility.mbmm(1.0/Je[...,:1].unsqueeze(-1), J[...,0,1:].unsqueeze(-2)) + erate = solvers.scalar_newton(RJ, torch.sign(csr) * self.guess_erate) + + ydot, J, Je, _ = self.model( + t, torch.cat([cs.unsqueeze(-1), y[..., 1:]], dim=-1), erate, cT + ) + + # There is an annoying extra term that is the derivative of the history rate with respect to the + # solved for strain rate times the derivative of the strain rate with respect to history + t1 = Je[..., 1:].unsqueeze(-1) + t2 = utility.mbmm(1.0 / Je[..., :1].unsqueeze(-1), J[..., 0, 1:].unsqueeze(-2)) t3 = utility.mbmm(t1, t2) # Corrected jacobian - row1 = torch.cat([torch.zeros_like(J[...,0,0]).unsqueeze(-1), -J[..., 0, 1:] / Je[...,0][...,None]], dim = -1).unsqueeze(-2) - rest = torch.cat([torch.zeros_like(J[...,1:,0]).unsqueeze(-1), J[...,1:,1:] - t3], dim = -1) - jac = torch.cat([row1,rest], dim = -2) + row1 = torch.cat( + [ + torch.zeros_like(J[..., 0, 0]).unsqueeze(-1), + -J[..., 0, 1:] / Je[..., 0][..., None], + ], + dim=-1, + ).unsqueeze(-2) + rest = torch.cat( + [torch.zeros_like(J[..., 1:, 0]).unsqueeze(-1), J[..., 1:, 1:] - t3], dim=-1 + ) + jac = torch.cat([row1, rest], dim=-2) - return torch.cat([erate.unsqueeze(-1), ydot[...,1:]], dim = -1), jac + return torch.cat([erate.unsqueeze(-1), ydot[..., 1:]], dim=-1), jac diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index ba8940b..ce2398b 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -275,7 +275,7 @@ class FixedGridBlockSolver: rtol (float): relative tolerance for Newton's method atol (float): absolute tolerance for Newton's method miter (int): maximum number of Newton iterations - linesearch (bool): whether to use + linesearch (bool): whether to use linear_solve_method (str): method to solve batched sparse Ax = b, options are currently "direct" or "dense" direct_solve_method (str): method to use for the direct solver, options are @@ -297,7 +297,7 @@ def __init__( rtol=1.0e-6, atol=1.0e-8, miter=200, - linesearch = False, + linesearch=False, linear_solve_method="direct", direct_solve_method="thomas", direct_solve_min_size=0, @@ -305,7 +305,7 @@ def __init__( guess_type="zero", guess_history=None, throw_on_fail=False, - offset_step = 0, + offset_step=0, **kwargs, ): # Store basic info about the system @@ -383,14 +383,13 @@ def integrate(self, t, cache_adjoint=False): result[0] = self.y0 incs = self._gen_increments(t) - - for k1,k2 in zip(incs[:-1], incs[1:]): - result[k1 : k2] = self.block_update( - t[k1 : k2], + for k1, k2 in zip(incs[:-1], incs[1:]): + result[k1:k2] = self.block_update( + t[k1:k2], t[k1 - 1], result[k1 - 1], self.func, - self._initial_guess(result, k1, k2-k1), + self._initial_guess(result, k1, k2 - k1), ) # Store for the backward pass, if we're going to do that @@ -411,7 +410,7 @@ def _gen_increments(self, t): ntotal = t.shape[0] if self.offset_step > 0: steps += [self.offset_step + 1] - steps += list(range(steps[-1],ntotal, self.n))[1:] + [ntotal] + steps += list(range(steps[-1], ntotal, self.n))[1:] + [ntotal] return steps @@ -470,15 +469,13 @@ def rewind(self, output_grad): # Calculate starts at last gradient prev_adjoint = output_grad[0] - + # Generate reverse increments incs = [1 + self.t.shape[0] - r for r in self._gen_increments(self.t)[::-1]] for k1, k2 in zip(incs[:-1], incs[1:]): # Could also cache these of course - _, J = self.func( - self.t[k1 - 1 : k2], self.result[k1 - 1 : k2] - ) + _, J = self.func(self.t[k1 - 1 : k2], self.result[k1 - 1 : k2]) full_adjoint = self.scheme.update_adjoint( self.t[k1 - 1 : k2].diff(dim=0), @@ -520,7 +517,7 @@ def block_update(self, t, t_start, y_start, func, y_guess): """ # Various useful sizes n = t.shape[0] # Number of time steps to do at once - + def RJ(dy): # Make things into a more rational shape dy = dy.reshape(self.batch_size, n, self.prob_size).transpose(0, 1) @@ -546,7 +543,7 @@ def RJ(dy): atol=self.atol, miter=self.miter, throw_on_fail=self.throw_on_fail, - linesearch = self.linesearch + linesearch=self.linesearch, ) return dy.reshape(self.batch_size, n, self.prob_size).transpose( diff --git a/pyoptmat/optimize.py b/pyoptmat/optimize.py index 26a88bc..06eb0b0 100644 --- a/pyoptmat/optimize.py +++ b/pyoptmat/optimize.py @@ -124,7 +124,7 @@ class DeterministicModel(Module): ics (list(torch.tensor)): initial conditions to use for each parameter """ - def __init__(self, maker, names, ics, use_cached_guess = False): + def __init__(self, maker, names, ics, use_cached_guess=False): super().__init__() self.maker = maker @@ -159,7 +159,7 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control): exp_control (torch.tensor): stress/strain control flag """ if self.use_cached_guess: - model = self.maker(*self.get_params(), guess_history = self.cached_solution) + model = self.maker(*self.get_params(), guess_history=self.cached_solution) else: model = self.maker(*self.get_params()) @@ -198,7 +198,9 @@ class StatisticalModel(PyroModule): entry i represents the noise in test type i """ - def __init__(self, maker, names, locs, scales, eps, nan_num=False, use_cached_guess = False): + def __init__( + self, maker, names, locs, scales, eps, nan_num=False, use_cached_guess=False + ): super().__init__() self.maker = maker @@ -242,7 +244,7 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control, exp_results=None exp_results (torch.tensor): true results for conditioning """ if self.use_cached_guess: - model = self.maker(*self.get_params(), guess_history = self.cached_solution) + model = self.maker(*self.get_params(), guess_history=self.cached_solution) else: model = self.maker(*self.get_params()) @@ -337,7 +339,7 @@ def __init__( param_suffix="_param", include_noise=False, weights=None, - use_cached_guess = False + use_cached_guess=False, ): super().__init__() @@ -556,8 +558,9 @@ def forward(self, exp_data, exp_cycles, exp_types, exp_control, exp_results=None # Sample the bottom level parameters if self.use_cached_guess: bmodel = self.maker( - *self.sample_bot(), extra_params=self.get_extra_params(), - guess_history = self.cached_solution + *self.sample_bot(), + extra_params=self.get_extra_params(), + guess_history=self.cached_solution, ) else: bmodel = self.maker( diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index 8deed4a..8082e2c 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -7,13 +7,14 @@ import warnings import torch -def scalar_bisection(fn, a, b, atol = 1.0e-6, miter = 100): + +def scalar_bisection(fn, a, b, atol=1.0e-6, miter=100): """ Solve logically scalar equations with bisection Args: fn (function): function returning scalar residual and jacobian - a (torch.tensor): lower bound + a (torch.tensor): lower bound b (torch.tensor): upper bound Keyword Args: @@ -22,11 +23,11 @@ def scalar_bisection(fn, a, b, atol = 1.0e-6, miter = 100): """ Ra, _ = fn(a) Rb, _ = fn(b) - + if not torch.all((torch.sign(Ra) + torch.sign(Rb)) == 0): raise RuntimeError("Initial values do not bisect in bisection solver") - - c = (a+b) / 2.0 + + c = (a + b) / 2.0 Rc, _ = fn(c) for i in range(miter): @@ -37,13 +38,14 @@ def scalar_bisection(fn, a, b, atol = 1.0e-6, miter = 100): bc = torch.sign(Rb) == torch.sign(Rc) a[ac] = c[ac] b[bc] = c[bc] - - c = (a+b) / 2.0 - Rc, _ = fn((a+b) / 2.0) + + c = (a + b) / 2.0 + Rc, _ = fn((a + b) / 2.0) return c -def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100): + +def scalar_newton(fn, x0, atol=1.0e-6, miter=100): """ Solve logically scalar equations with Newton's method @@ -61,22 +63,25 @@ def scalar_newton(fn, x0, atol = 1.0e-6, miter = 100): for i in range(miter): if torch.all(torch.abs(R) < atol): break - + x = x - R / J - + R, J = fn(x) else: - warnings.warn("Scalar implicit solve did not succeed. Results may be inaccurate...") - + warnings.warn( + "Scalar implicit solve did not succeed. Results may be inaccurate..." + ) + return x -def scalar_bisection_newton(fn, a, b, atol = 1.0e-6, miter = 100, biter = 10): + +def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10): """ Solve logically scalar equations by switching from bisection to Newton's method Args: fn (function): function returning scalar residual and jacobian - a (torch.tensor): lower bound + a (torch.tensor): lower bound b (torch.tensor): upper bound Keyword Args: @@ -84,8 +89,9 @@ def scalar_bisection_newton(fn, a, b, atol = 1.0e-6, miter = 100, biter = 10): biter: initial number of bisection iterations miter (int): max number of iterations for Newton's method """ - x = scalar_bisection(fn, a, b, atol = atol, miter = biter) - return scalar_newton(fn, x, atol = atol, miter = miter) + x = scalar_bisection(fn, a, b, atol=atol, miter=biter) + return scalar_newton(fn, x, atol=atol, miter=miter) + def newton_raphson_bt( fn, x0, linsolver="lu", rtol=1e-6, atol=1e-10, miter=100, max_bt=5 diff --git a/test/test_flowrules.py b/test/test_flowrules.py index 65c54dd..1b8f735 100644 --- a/test/test_flowrules.py +++ b/test/test_flowrules.py @@ -209,6 +209,7 @@ def setUp(self): self.rtol = 1e-4 self.rtol2 = 1e-4 + class TestPerfectRateIndependentPlasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch ): @@ -230,6 +231,7 @@ def setUp(self): self.rtol = 1e-4 self.rtol2 = 1e-4 + class TestWrappedRIIsoKinViscoplasticity( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch ): @@ -586,6 +588,7 @@ def fn(i): self.assertTrue(np.allclose(i1, i2, rtol=1.0e-4)) + class TestSuperimposedFlowRate( unittest.TestCase, CommonFlowRule, CommonFlowRuleBatchBatch ): @@ -734,7 +737,7 @@ def setUp(self): self.kin = hardening.ChabocheHardeningModel(CP(self.C), CP(self.g)) self.model = flowrules.IsoKinRateIndependentPlasticity( - CP(self.E), CP(self.sy), self.iso, self.kin, s = 0.01 + CP(self.E), CP(self.sy), self.iso, self.kin, s=0.01 ) self.s = torch.linspace(160, 240, self.nbatch) @@ -763,52 +766,88 @@ def setUp(self): self.ep_guess = self.erate.clone() def test_d_residual_dep(self): - exact = self.model.ep_jacobian_ep(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) - numer = utility.batch_differentiate(lambda x: self.model.ep_residual(x, self.s, self.h, self.t, self.T, self.erate), self.ep_guess).unsqueeze(-1) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + exact = self.model.ep_jacobian_ep( + self.ep_guess, self.s, self.h, self.t, self.T, self.erate + ) + numer = utility.batch_differentiate( + lambda x: self.model.ep_residual( + x, self.s, self.h, self.t, self.T, self.erate + ), + self.ep_guess, + ).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_d_residual_ds(self): - exact = self.model.ep_jacobian_s(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) - numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, x, self.h, self.t, self.T, self.erate), self.s).unsqueeze(-1) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + exact = self.model.ep_jacobian_s( + self.ep_guess, self.s, self.h, self.t, self.T, self.erate + ) + numer = utility.batch_differentiate( + lambda x: self.model.ep_residual( + self.ep_guess, x, self.h, self.t, self.T, self.erate + ), + self.s, + ).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_d_residual_dh(self): - exact = self.model.ep_jacobian_h(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) - numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, x, self.t, self.T, self.erate), self.h) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + exact = self.model.ep_jacobian_h( + self.ep_guess, self.s, self.h, self.t, self.T, self.erate + ) + numer = utility.batch_differentiate( + lambda x: self.model.ep_residual( + self.ep_guess, self.s, x, self.t, self.T, self.erate + ), + self.h, + ) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_d_residual_de(self): - exact = self.model.ep_jacobian_e(self.ep_guess, self.s, self.h, self.t, self.T, self.erate) - numer = utility.batch_differentiate(lambda x: self.model.ep_residual(self.ep_guess, self.s, self.h, self.t, self.T, x), self.erate).unsqueeze(-1) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + exact = self.model.ep_jacobian_e( + self.ep_guess, self.s, self.h, self.t, self.T, self.erate + ) + numer = utility.batch_differentiate( + lambda x: self.model.ep_residual( + self.ep_guess, self.s, self.h, self.t, self.T, x + ), + self.erate, + ).unsqueeze(-1) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_d_sig(self): x = torch.linspace(-1, 1, 10) exact = self.model.dsig(x) numer = utility.batch_differentiate(lambda x: self.model.sig(x), x) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_df_ds(self): exact = self.model.df_ds(self.s, self.h, self.T) - numer = utility.batch_differentiate(lambda x: self.model.f(x, self.h, self.T), self.s) + numer = utility.batch_differentiate( + lambda x: self.model.f(x, self.h, self.T), self.s + ) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_dmix_fn_ds(self): exact = self.model.dmix_fn_ds(self.s, self.h, self.T) - numer = utility.batch_differentiate(lambda x: self.model.mix_fn(x, self.h, self.T), self.s) - - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol, atol = 1.0e-4)) + numer = utility.batch_differentiate( + lambda x: self.model.mix_fn(x, self.h, self.T), self.s + ) + + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol, atol=1.0e-4)) def test_df_dh(self): exact = self.model.df_dh(self.s, self.h, self.T) - numer = utility.batch_differentiate(lambda x: self.model.f(self.s, x, self.T), self.h) + numer = utility.batch_differentiate( + lambda x: self.model.f(self.s, x, self.T), self.h + ) - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) def test_dmix_fn_dh(self): exact = self.model.dmix_fn_dh(self.s, self.h, self.T) - numer = utility.batch_differentiate(lambda x: self.model.mix_fn(self.s, x, self.T), self.h) - - self.assertTrue(torch.allclose(exact, numer, rtol = self.rtol)) + numer = utility.batch_differentiate( + lambda x: self.model.mix_fn(self.s, x, self.T), self.h + ) + + self.assertTrue(torch.allclose(exact, numer, rtol=self.rtol)) diff --git a/test/test_models.py b/test/test_models.py index 8db3066..ec083c9 100644 --- a/test/test_models.py +++ b/test/test_models.py @@ -71,7 +71,7 @@ def test_derivs_stress(self): ddv = utility.batch_differentiate( lambda x: use.forward(self.t, x)[0], self.state_stress ) - + self.assertTrue(np.allclose(dv, ddv, rtol=1e-4, atol=1e-4)) def test_partial_state(self): @@ -327,6 +327,7 @@ def setUp(self): self.t = self.times[2] self.step = 2 + class TestIsoKinRIPlasticity(unittest.TestCase, CommonModel, CommonModelBatchBatch): def setUp(self): self.E = torch.tensor(100000.0) @@ -362,8 +363,8 @@ def setUp(self): self.state_strain = ( torch.tensor([[90.0, 30.0, 10.0], [100.0, 10.0, 15.0], [101.0, 50.0, 60.0]]) ) / 3.0 - self.state_stress = ( - torch.tensor([[0.05, 30.0, 10.0], [0.07, 10.0, 15.0], [0.08, 50.0, 60.0]]) + self.state_stress = torch.tensor( + [[0.05, 30.0, 10.0], [0.07, 10.0, 15.0], [0.08, 50.0, 60.0]] ) self.t = self.times[2] @@ -534,4 +535,3 @@ def setUp(self): self.t = self.times[2] self.step = 2 - diff --git a/test/test_ode.py b/test/test_ode.py index 399ebb1..267c28f 100644 --- a/test/test_ode.py +++ b/test/test_ode.py @@ -181,6 +181,7 @@ def test_backward(self): self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) + class TestFallingBatchUnequal(unittest.TestCase): def setUp(self): self.m = 1.0 / 4.0 @@ -206,7 +207,7 @@ def test_forward(self): self.times, method="forward-euler", block_size=self.nchunk, - offset_step = 2 + offset_step=2, ) self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) @@ -219,7 +220,7 @@ def test_backward(self): self.times, method="backward-euler", block_size=self.nchunk, - offset_step = 2 + offset_step=2, ) self.assertTrue(torch.max(torch.abs((numerical[1:] - exact[1:]))) < 1.0e-1) From c50b0ec24e32b007fcc32935ad3bf914c10dca2a Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 11 Jan 2024 21:28:28 -0600 Subject: [PATCH 12/18] Started linting --- pyoptmat/chunktime.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyoptmat/chunktime.py b/pyoptmat/chunktime.py index 43150d1..3836f95 100644 --- a/pyoptmat/chunktime.py +++ b/pyoptmat/chunktime.py @@ -58,10 +58,9 @@ def newton_raphson_chunk( while (i < miter) and torch.any( torch.logical_not(torch.logical_or(nR <= atol, nR / nR0 <= rtol)) ): - print(i, torch.max(nR)) dx = solver.solve(J, R) if linesearch: - x, R, J, nR, alpha = chunk_linesearch(x, dx, fn, R, rtol, atol) + x, R, J, nR = chunk_linesearch(x, dx, fn, R, rtol, atol) else: x -= dx R, J = fn(x) @@ -82,7 +81,8 @@ def chunk_linesearch( """ Backtracking linesearch for the chunk NR algorithm. - Terminates when the Armijo criteria is reached, or you exceed some maximum iterations, or when you would meet the + Terminates when the Armijo criteria is reached, or you exceed + some maximum iterations, or when you would meet the convergence requirements with the current alpha Args: @@ -114,7 +114,7 @@ def chunk_linesearch( break alpha[nR >= crit] /= sigma - return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim=-1), alpha + return x - dx * alpha.unsqueeze(-1), R, J, torch.norm(R, dim=-1) class BidiagonalOperator(torch.nn.Module): From e06af04720fe010347c56e9562419ffaa1246207 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Fri, 12 Jan 2024 19:48:57 -0600 Subject: [PATCH 13/18] Okay, almost done --- pyoptmat/models.py | 14 ++++++++++++-- pyoptmat/solvers.py | 8 +++++--- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 13a0f93..c475741 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -283,10 +283,11 @@ class ModelIntegrator(nn.Module): """ - def __init__(self, model, *args, use_adjoint=True, bisect_first=False, **kwargs): + def __init__(self, model, *args, use_adjoint=True, bisect_first=False, throw_on_scalar_fail = False, **kwargs): super().__init__(*args) self.model = model self.use_adjoint = use_adjoint + self.throw_on_scalar_fail = throw_on_scalar_fail self.kwargs_for_integration = kwargs if self.use_adjoint: @@ -325,6 +326,7 @@ def solve_both(self, times, temperatures, idata, control): temperatures, control, bisect_first=self.bisect_first, + throw_on_scalar_fail = self.throw_on_scalar_fail ) return self.imethod(bmodel, init, times, **self.kwargs_for_integration) @@ -405,6 +407,8 @@ def solve_stress(self, times, stresses, temperatures): stress_rate_interpolator, stress_interpolator, temperature_interpolator, + bisect_first = self.bisect_first, + throw_on_scalar_fail = self.throw_on_scalar_fail ) return self.imethod(smodel, init, times, **self.kwargs_for_integration) @@ -442,6 +446,7 @@ def __init__( temps, control, bisect_first=False, + throw_on_scalar_fail = False, *args, **kwargs ): @@ -473,6 +478,7 @@ def __init__( times[..., self.scontrol], temps[..., self.scontrol] ), bisect_first=bisect_first, + throw_on_scalar_fail = throw_on_scalar_fail ) def forward(self, t, y): @@ -556,6 +562,7 @@ def __init__( max_erate=1e3, guess_erate=1.0e-3, bisect_first=False, + throw_on_scalar_fail = False, *args, **kwargs ): @@ -568,6 +575,7 @@ def __init__( self.max_erate = max_erate self.bisect_first = bisect_first self.guess_erate = guess_erate + self.throw_on_scalar_fail = throw_on_scalar_fail def forward(self, t, y): """ @@ -597,9 +605,11 @@ def RJ(erate): RJ, torch.ones_like(y[..., 0]) * self.min_erate, torch.ones_like(y[..., 0]) * self.max_erate, + throw_on_fail = self.throw_on_scalar_fail ) else: - erate = solvers.scalar_newton(RJ, torch.sign(csr) * self.guess_erate) + erate = solvers.scalar_newton(RJ, torch.sign(csr) * self.guess_erate, + throw_on_fail = self.throw_on_scalar_fail) ydot, J, Je, _ = self.model( t, torch.cat([cs.unsqueeze(-1), y[..., 1:]], dim=-1), erate, cT diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index 8082e2c..96a9600 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -45,7 +45,7 @@ def scalar_bisection(fn, a, b, atol=1.0e-6, miter=100): return c -def scalar_newton(fn, x0, atol=1.0e-6, miter=100): +def scalar_newton(fn, x0, atol=1.0e-6, miter=100, throw_on_fail = False): """ Solve logically scalar equations with Newton's method @@ -68,6 +68,8 @@ def scalar_newton(fn, x0, atol=1.0e-6, miter=100): R, J = fn(x) else: + if throw_on_fail: + raise RuntimeError("Scalar solve failed") warnings.warn( "Scalar implicit solve did not succeed. Results may be inaccurate..." ) @@ -75,7 +77,7 @@ def scalar_newton(fn, x0, atol=1.0e-6, miter=100): return x -def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10): +def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10, throw_on_fail = False): """ Solve logically scalar equations by switching from bisection to Newton's method @@ -90,7 +92,7 @@ def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10): miter (int): max number of iterations for Newton's method """ x = scalar_bisection(fn, a, b, atol=atol, miter=biter) - return scalar_newton(fn, x, atol=atol, miter=miter) + return scalar_newton(fn, x, atol=atol, miter=miter, throw_on_fail = throw_on_fail) def newton_raphson_bt( From b2027e990f1b690e3c21d2caf0c9949a2d7e4cac Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Tue, 16 Jan 2024 22:14:51 -0600 Subject: [PATCH 14/18] KM model should be scaled --- pyoptmat/flowrules.py | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index 6073623..306b090 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -402,7 +402,7 @@ class SoftKocksMeckingRegimeFlowRule(FlowRule): Args: model1 (flowrules.FlowRule): first flow rule model2 (flowrules.FlowRule): second flow rule - g0 (torch.tensor): activation energy threshold + A, B, C (torch.tensor): Kocks-Mecking parameters mu (temperature.TemperatureParameter): shear modulus b (torch.tensor): burgers vector eps0 (torch.tensor): reference strain rate @@ -410,30 +410,35 @@ class SoftKocksMeckingRegimeFlowRule(FlowRule): sf (torch.tensor): sharpness parameter Keyword Args: + A_scale, B_scale, C_scale: scaling functions for K-M parameters eps (float): default 1e-20, offset to avoid divide-by-zero - g0_scale (function): scaling function for g0, - defaults to no scaling """ def __init__( self, model1, model2, - g0, + A, + B, + C, mu, b, eps0, k, sf, + A_scale = lambda x: x, + B_scale = lambda x: x, + C_scale = lambda x: x, eps=torch.tensor(1e-20), - g0_scale=lambda x: x, ): super().__init__() self.model1 = model1 self.model2 = model2 - self.g0 = g0 + self.A = A + self.B = B + self.C = C self.mu = mu self.b = b @@ -443,7 +448,9 @@ def __init__( self.sf = sf self.eps = eps - self.g0_scale = g0_scale + self.A_scale = A_scale + self.B_scale = B_scale + self.C_scale = C_scale # Check for conformal history vectors if self.model1.nhist != self.model2.nhist: @@ -505,9 +512,12 @@ def f(self, T, e): torch.tensor: value of the weighting function """ return ( - torch.tanh(self.sf * (self.g(T, e) - self.g0_scale(self.g0))) + 1.0 + torch.tanh(self.sf * (self.g(T, e) - self.g0())) + 1.0 ) / 2.0 + def g0(self): + return (self.C_scale.scale(self.C) - self.B_scale.scale(self.B)) / self.A_scale(self.A) + def df_e(self, T, e): """ The derivative of the weight function with respect to @@ -524,7 +534,7 @@ def df_e(self, T, e): self.sf / ( 2.0 - * torch.cosh(self.sf * (self.g(T, e) - self.g0_scale(self.g0))) ** 2.0 + * torch.cosh(self.sf * (self.g(T, e) - self.g0())) ** 2.0 ) * self.dg_e(T, e) ) From 45acf99ff78873ca06958ef5b30ee7a819be79e0 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 21 Mar 2024 16:57:58 -0500 Subject: [PATCH 15/18] Blacking --- pyoptmat/flowrules.py | 19 ++++++++----------- pyoptmat/models.py | 31 +++++++++++++++++++++---------- pyoptmat/solvers.py | 8 +++++--- 3 files changed, 34 insertions(+), 24 deletions(-) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index 306b090..987e2dd 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -427,9 +427,9 @@ def __init__( eps0, k, sf, - A_scale = lambda x: x, - B_scale = lambda x: x, - C_scale = lambda x: x, + A_scale=lambda x: x, + B_scale=lambda x: x, + C_scale=lambda x: x, eps=torch.tensor(1e-20), ): super().__init__() @@ -511,12 +511,12 @@ def f(self, T, e): Returns: torch.tensor: value of the weighting function """ - return ( - torch.tanh(self.sf * (self.g(T, e) - self.g0())) + 1.0 - ) / 2.0 + return (torch.tanh(self.sf * (self.g(T, e) - self.g0())) + 1.0) / 2.0 def g0(self): - return (self.C_scale.scale(self.C) - self.B_scale.scale(self.B)) / self.A_scale(self.A) + return (self.C_scale.scale(self.C) - self.B_scale.scale(self.B)) / self.A_scale( + self.A + ) def df_e(self, T, e): """ @@ -532,10 +532,7 @@ def df_e(self, T, e): """ return ( self.sf - / ( - 2.0 - * torch.cosh(self.sf * (self.g(T, e) - self.g0())) ** 2.0 - ) + / (2.0 * torch.cosh(self.sf * (self.g(T, e) - self.g0())) ** 2.0) * self.dg_e(T, e) ) diff --git a/pyoptmat/models.py b/pyoptmat/models.py index c475741..1aa9e64 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -283,7 +283,15 @@ class ModelIntegrator(nn.Module): """ - def __init__(self, model, *args, use_adjoint=True, bisect_first=False, throw_on_scalar_fail = False, **kwargs): + def __init__( + self, + model, + *args, + use_adjoint=True, + bisect_first=False, + throw_on_scalar_fail=False, + **kwargs + ): super().__init__(*args) self.model = model self.use_adjoint = use_adjoint @@ -326,7 +334,7 @@ def solve_both(self, times, temperatures, idata, control): temperatures, control, bisect_first=self.bisect_first, - throw_on_scalar_fail = self.throw_on_scalar_fail + throw_on_scalar_fail=self.throw_on_scalar_fail, ) return self.imethod(bmodel, init, times, **self.kwargs_for_integration) @@ -407,8 +415,8 @@ def solve_stress(self, times, stresses, temperatures): stress_rate_interpolator, stress_interpolator, temperature_interpolator, - bisect_first = self.bisect_first, - throw_on_scalar_fail = self.throw_on_scalar_fail + bisect_first=self.bisect_first, + throw_on_scalar_fail=self.throw_on_scalar_fail, ) return self.imethod(smodel, init, times, **self.kwargs_for_integration) @@ -446,7 +454,7 @@ def __init__( temps, control, bisect_first=False, - throw_on_scalar_fail = False, + throw_on_scalar_fail=False, *args, **kwargs ): @@ -478,7 +486,7 @@ def __init__( times[..., self.scontrol], temps[..., self.scontrol] ), bisect_first=bisect_first, - throw_on_scalar_fail = throw_on_scalar_fail + throw_on_scalar_fail=throw_on_scalar_fail, ) def forward(self, t, y): @@ -562,7 +570,7 @@ def __init__( max_erate=1e3, guess_erate=1.0e-3, bisect_first=False, - throw_on_scalar_fail = False, + throw_on_scalar_fail=False, *args, **kwargs ): @@ -605,11 +613,14 @@ def RJ(erate): RJ, torch.ones_like(y[..., 0]) * self.min_erate, torch.ones_like(y[..., 0]) * self.max_erate, - throw_on_fail = self.throw_on_scalar_fail + throw_on_fail=self.throw_on_scalar_fail, ) else: - erate = solvers.scalar_newton(RJ, torch.sign(csr) * self.guess_erate, - throw_on_fail = self.throw_on_scalar_fail) + erate = solvers.scalar_newton( + RJ, + torch.sign(csr) * self.guess_erate, + throw_on_fail=self.throw_on_scalar_fail, + ) ydot, J, Je, _ = self.model( t, torch.cat([cs.unsqueeze(-1), y[..., 1:]], dim=-1), erate, cT diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index 96a9600..7c7a1c1 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -45,7 +45,7 @@ def scalar_bisection(fn, a, b, atol=1.0e-6, miter=100): return c -def scalar_newton(fn, x0, atol=1.0e-6, miter=100, throw_on_fail = False): +def scalar_newton(fn, x0, atol=1.0e-6, miter=100, throw_on_fail=False): """ Solve logically scalar equations with Newton's method @@ -77,7 +77,9 @@ def scalar_newton(fn, x0, atol=1.0e-6, miter=100, throw_on_fail = False): return x -def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10, throw_on_fail = False): +def scalar_bisection_newton( + fn, a, b, atol=1.0e-6, miter=100, biter=10, throw_on_fail=False +): """ Solve logically scalar equations by switching from bisection to Newton's method @@ -92,7 +94,7 @@ def scalar_bisection_newton(fn, a, b, atol=1.0e-6, miter=100, biter=10, throw_on miter (int): max number of iterations for Newton's method """ x = scalar_bisection(fn, a, b, atol=atol, miter=biter) - return scalar_newton(fn, x, atol=atol, miter=miter, throw_on_fail = throw_on_fail) + return scalar_newton(fn, x, atol=atol, miter=miter, throw_on_fail=throw_on_fail) def newton_raphson_bt( From 4f611f56947c3c5ffc629ed7e5085c08be3af0f4 Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 21 Mar 2024 20:33:45 -0500 Subject: [PATCH 16/18] Hopefully formatting fixed --- pyoptmat/flowrules.py | 5 ++++- pyoptmat/models.py | 9 +++++---- pyoptmat/ode.py | 2 +- pyoptmat/optimize.py | 2 ++ pyoptmat/solvers.py | 4 ++-- 5 files changed, 14 insertions(+), 8 deletions(-) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index 987e2dd..ee05536 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -1,4 +1,4 @@ -# pylint: disable=abstract-method, useless-super-delegation, line-too-long, duplicate-code, too-many-lines +# pylint: disable=abstract-method, useless-super-delegation, line-too-long, duplicate-code, too-many-lines, too-many-arguments, too-many-public-methods """ Module containing inelastic flow rules. These provide the rate of the @@ -514,6 +514,9 @@ def f(self, T, e): return (torch.tanh(self.sf * (self.g(T, e) - self.g0())) + 1.0) / 2.0 def g0(self): + """ + The intercept value + """ return (self.C_scale.scale(self.C) - self.B_scale.scale(self.B)) / self.A_scale( self.A ) diff --git a/pyoptmat/models.py b/pyoptmat/models.py index 1aa9e64..ccca4dc 100644 --- a/pyoptmat/models.py +++ b/pyoptmat/models.py @@ -453,9 +453,9 @@ def __init__( base, temps, control, + *args, bisect_first=False, throw_on_scalar_fail=False, - *args, **kwargs ): super().__init__(*args, **kwargs) @@ -566,12 +566,12 @@ def __init__( srate_fn, stress_fn, T_fn, + *args, min_erate=-1e2, max_erate=1e3, guess_erate=1.0e-3, bisect_first=False, throw_on_scalar_fail=False, - *args, **kwargs ): super().__init__(*args, **kwargs) @@ -626,8 +626,9 @@ def RJ(erate): t, torch.cat([cs.unsqueeze(-1), y[..., 1:]], dim=-1), erate, cT ) - # There is an annoying extra term that is the derivative of the history rate with respect to the - # solved for strain rate times the derivative of the strain rate with respect to history + # There is an annoying extra term that is the derivative of the + # history rate with respect to the solved for strain rate times + # the derivative of the strain rate with respect to history t1 = Je[..., 1:].unsqueeze(-1) t2 = utility.mbmm(1.0 / Je[..., :1].unsqueeze(-1), J[..., 0, 1:].unsqueeze(-2)) t3 = utility.mbmm(t1, t2) diff --git a/pyoptmat/ode.py b/pyoptmat/ode.py index ce2398b..76ef0d8 100644 --- a/pyoptmat/ode.py +++ b/pyoptmat/ode.py @@ -1,4 +1,4 @@ -# pylint: disable=too-many-arguments +# pylint: disable=too-many-arguments, too-many-instance-attributes """ Module defining the key objects and functions to integrate ODEs diff --git a/pyoptmat/optimize.py b/pyoptmat/optimize.py index 06eb0b0..5bd6dcb 100644 --- a/pyoptmat/optimize.py +++ b/pyoptmat/optimize.py @@ -1,3 +1,5 @@ +# pylint: disable=too-many-instance-attributes, too-many-arguments, not-context-manager + """ Objects and helper functions to help with deterministic model calibration and statistical inference. diff --git a/pyoptmat/solvers.py b/pyoptmat/solvers.py index 7c7a1c1..658af48 100644 --- a/pyoptmat/solvers.py +++ b/pyoptmat/solvers.py @@ -30,7 +30,7 @@ def scalar_bisection(fn, a, b, atol=1.0e-6, miter=100): c = (a + b) / 2.0 Rc, _ = fn(c) - for i in range(miter): + for _ in range(miter): if torch.all(torch.abs(Rc) < atol): break @@ -60,7 +60,7 @@ def scalar_newton(fn, x0, atol=1.0e-6, miter=100, throw_on_fail=False): x = x0 R, J = fn(x) - for i in range(miter): + for _ in range(miter): if torch.all(torch.abs(R) < atol): break From 65fcd414ee99a38714be031ad141d57c48ba8bee Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 21 Mar 2024 20:49:23 -0500 Subject: [PATCH 17/18] Fixed tests --- pyoptmat/flowrules.py | 2 +- test/test_flowrules.py | 18 ++++++++++++++---- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index ee05536..cd873c9 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -517,7 +517,7 @@ def g0(self): """ The intercept value """ - return (self.C_scale.scale(self.C) - self.B_scale.scale(self.B)) / self.A_scale( + return (self.C_scale(self.C) - self.B_scale(self.B)) / self.A_scale( self.A ) diff --git a/test/test_flowrules.py b/test/test_flowrules.py index 1b8f735..ed8cc49 100644 --- a/test/test_flowrules.py +++ b/test/test_flowrules.py @@ -407,13 +407,18 @@ def setUp(self): self.eps0 = 1.0 self.k = 1.0 - self.g0 = 1.5 + self.A = torch.tensor(-3.35) + self.B = torch.tensor(-3.23) + self.C = torch.tensor(-5.82) + self.sf = 10.0 self.model = flowrules.SoftKocksMeckingRegimeFlowRule( self.model1, self.model2, - self.g0, + self.A, + self.B, + self.C, self.mu, self.b, self.eps0, @@ -486,13 +491,18 @@ def setUp(self): self.eps0 = 1.0 self.k = 1.0 - self.g0 = 1.5 + self.A = torch.tensor(-3.35) + self.B = torch.tensor(-3.23) + self.C = torch.tensor(-5.82) + self.sf = 10.0 self.model = flowrules.SoftKocksMeckingRegimeFlowRule( self.model1, self.model2, - self.g0, + self.A, + self.B, + self.C, self.mu, self.b, self.eps0, From 8aedafb463a551a4cb12fcbbeb89b4a46b075a8f Mon Sep 17 00:00:00 2001 From: Mark Messner Date: Thu, 21 Mar 2024 21:03:18 -0500 Subject: [PATCH 18/18] Fix examples --- examples/ode/damping.py | 7 ++----- .../kocks-mecking-differentiable.py | 3 +-- pyoptmat/flowrules.py | 4 +--- 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/examples/ode/damping.py b/examples/ode/damping.py index 61bcc7f..e714dbe 100755 --- a/examples/ode/damping.py +++ b/examples/ode/damping.py @@ -35,14 +35,11 @@ sys.path.append("../..") from pyoptmat import ode, utility, neuralode -if torch.cuda.is_available(): - dev = "cuda:0" -else: - dev = "cpu" +dev = "cpu" device = torch.device(dev) # I'm assuming single precision will be fine for this but for now use doubles -torch.set_default_tensor_type(torch.DoubleTensor) +torch.set_default_dtype(torch.float64) class MassDamperSpring(torch.nn.Module): diff --git a/examples/structural-material-models/kocks-mecking-differentiable.py b/examples/structural-material-models/kocks-mecking-differentiable.py index b0d15f4..1712b3f 100755 --- a/examples/structural-material-models/kocks-mecking-differentiable.py +++ b/examples/structural-material-models/kocks-mecking-differentiable.py @@ -48,7 +48,6 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): mu = temperature.PolynomialScaling( [-2.60610204e-05, 3.61911162e-02, -4.23765368e01, 8.44212545e04] ) - g0 = torch.tensor(0.771) k = torch.tensor(1.38064e-20) b = torch.tensor(2.019e-7) eps0 = torch.tensor(1.0e6) @@ -81,7 +80,7 @@ def calculate_yield(strain, stress, offset=0.2 / 100.0): sf = torch.tensor(50.0) flowrule = flowrules.SoftKocksMeckingRegimeFlowRule( - ri_flowrule, rd_flowrule, g0, mu, b, eps0, k, sf + ri_flowrule, rd_flowrule, A, B, C, mu, b, eps0, k, sf ) model = models.InelasticModel(E, flowrule) diff --git a/pyoptmat/flowrules.py b/pyoptmat/flowrules.py index cd873c9..a422325 100644 --- a/pyoptmat/flowrules.py +++ b/pyoptmat/flowrules.py @@ -517,9 +517,7 @@ def g0(self): """ The intercept value """ - return (self.C_scale(self.C) - self.B_scale(self.B)) / self.A_scale( - self.A - ) + return (self.C_scale(self.C) - self.B_scale(self.B)) / self.A_scale(self.A) def df_e(self, T, e): """