diff --git a/examples/2_Intermediate/strain_optimization.py b/examples/2_Intermediate/strain_optimization.py index 5087d4fd4..cb04ba2b7 100755 --- a/examples/2_Intermediate/strain_optimization.py +++ b/examples/2_Intermediate/strain_optimization.py @@ -13,7 +13,7 @@ import numpy as np from scipy.optimize import minimize -from simsopt.geo import StrainOpt, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty +from simsopt.geo import CoilStrain, LPTorsionalStrainPenalty, LPBinormalCurvatureStrainPenalty from simsopt.geo import FrameRotation, FramedCurveFrenet, CurveXYZFourier from simsopt.configs import get_hsx_data from simsopt.util import in_github_actions @@ -40,7 +40,7 @@ Jbin = LPBinormalCurvatureStrainPenalty( framedcurve, p=2, threshold=cur_threshold) -strain = StrainOpt(framedcurve, width) +strain = CoilStrain(framedcurve, width) JF = Jtor + Jbin diff --git a/src/simsopt/geo/__init__.py b/src/simsopt/geo/__init__.py index 98d1a8129..bc6d0efe2 100644 --- a/src/simsopt/geo/__init__.py +++ b/src/simsopt/geo/__init__.py @@ -22,7 +22,6 @@ from .surfacexyzfourier import * from .surfacexyztensorfourier import * from .strain_optimization import * -from .framedcurve import * from .permanent_magnet_grid import * diff --git a/src/simsopt/geo/finitebuild.py b/src/simsopt/geo/finitebuild.py index 3d1783e9d..87a4cc60b 100644 --- a/src/simsopt/geo/finitebuild.py +++ b/src/simsopt/geo/finitebuild.py @@ -1,12 +1,5 @@ import numpy as np -import jax.numpy as jnp -from jax import vjp, jvp - -import simsoptpp as sopp -from .._core.optimizable import Optimizable -from .._core.derivative import Derivative -from .curve import Curve -from .jit import jit + from .framedcurve import FramedCurve, FrameRotation, ZeroRotation, FramedCurveCentroid, FramedCurveFrenet """ @@ -50,7 +43,6 @@ def recompute_bell(self, parent=None): def gamma_impl(self, gamma, quadpoints): assert quadpoints.shape[0] == self.curve.quadpoints.shape[0] assert np.linalg.norm(quadpoints - self.curve.quadpoints) < 1e-15 - c = self.curve t, n, b = self.framedcurve.rotated_frame() gamma[:] = self.curve.gamma() + self.dn * n + self.db * b diff --git a/src/simsopt/geo/framedcurve.py b/src/simsopt/geo/framedcurve.py index d580d5e4b..dc883418f 100644 --- a/src/simsopt/geo/framedcurve.py +++ b/src/simsopt/geo/framedcurve.py @@ -261,7 +261,6 @@ def frame_torsion(self): gamma = self.curve.gamma() d1gamma = self.curve.gammadash() d2gamma = self.curve.gammadashdash() - d3gamma = self.curve.gammadashdashdash() alpha = self.rotation.alpha(self.curve.quadpoints) alphadash = self.rotation.alphadash(self.curve.quadpoints) return self.torsion(gamma, d1gamma, d2gamma, alpha, alphadash) @@ -270,7 +269,6 @@ def frame_binormal_curvature(self): gamma = self.curve.gamma() d1gamma = self.curve.gammadash() d2gamma = self.curve.gammadashdash() - d3gamma = self.curve.gammadashdashdash() alpha = self.rotation.alpha(self.curve.quadpoints) alphadash = self.rotation.alphadash(self.curve.quadpoints) return self.binorm(gamma, d1gamma, d2gamma, alpha, alphadash) @@ -422,7 +420,7 @@ def __init__(self, quadpoints): .. code-block:: python - rot = FilamentRotation(...) + rot = FrameRotation(...) rot.fix_all() """ @@ -604,12 +602,6 @@ def inner(a, b): return np.sum(a*b, axis=1) -torsion2vjp0 = jit(lambda ndash, b, v: vjp( - lambda nd: torsion_pure(nd, b), ndash)[1](v)[0]) -torsion2vjp1 = jit(lambda ndash, b, v: vjp( - lambda bi: torsion_pure(ndash, bi), b)[1](v)[0]) - - def torsion_pure_frenet(gamma, gammadash, gammadashdash, gammadashdashdash, alpha, alphadash): """Torsion function for export/evaulate coil sets""" diff --git a/src/simsopt/geo/strain_optimization.py b/src/simsopt/geo/strain_optimization.py index cbe7ac5d9..59400c107 100644 --- a/src/simsopt/geo/strain_optimization.py +++ b/src/simsopt/geo/strain_optimization.py @@ -6,7 +6,7 @@ from simsopt.geo.curveobjectives import Lp_torsion_pure __all__ = ['LPBinormalCurvatureStrainPenalty', - 'LPTorsionalStrainPenalty', 'StrainOpt'] + 'LPTorsionalStrainPenalty', 'CoilStrain'] class LPBinormalCurvatureStrainPenalty(Optimizable): @@ -15,14 +15,21 @@ class LPBinormalCurvatureStrainPenalty(Optimizable): of the binormal curvature strain, and penalizes where the local strain exceeds a threshold .. math:: - J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{bend}} - \epsilon_0, 0)^p ~dl + J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{bend}} - \epsilon_0, 0)^p ~dl, - where :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``. + where + + .. math:: + \epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2}, + + :math:`w` is the width of the tape, :math:`\hat{\textbf{b}}` is the + frame binormal vector, :math:`\boldsymbol{\kappa}` is the curvature vector of the + filamentary coil, and :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``. """ def __init__(self, framedcurve, width=1e-3, p=2, threshold=0): self.framedcurve = framedcurve - self.strain = StrainOpt(framedcurve, width) + self.strain = CoilStrain(framedcurve, width) self.width = width self.p = p self.threshold = threshold @@ -65,12 +72,18 @@ class LPTorsionalStrainPenalty(Optimizable): .. math:: J = \frac{1}{p} \int_{\text{curve}} \text{max}(\epsilon_{\text{tor}} - \epsilon_0, 0)^p ~dl - where :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``. + where + + .. math:: + \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12}, + + :math:`\tau` is the torsion of the tape frame, :math:`w` is the width of the tape, + and :math:`\epsilon_0` is a threshold strain, given by the argument ``threshold``. """ def __init__(self, framedcurve, width=1e-3, p=2, threshold=0): self.framedcurve = framedcurve - self.strain = StrainOpt(framedcurve, width) + self.strain = CoilStrain(framedcurve, width) self.width = width self.p = p self.threshold = threshold @@ -105,7 +118,7 @@ def dJ(self): return_fn_map = {'J': J, 'dJ': dJ} -class StrainOpt(Optimizable): +class CoilStrain(Optimizable): r""" This class evaluates the torsional and binormal curvature strains on HTS, based on a filamentary model of the coil and the orientation of the HTS tape. @@ -119,13 +132,18 @@ class StrainOpt(Optimizable): the expressions for the strains are: .. math:: - \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12} + \epsilon_{\text{tor}} = \frac{\tau^2 w^2}{12} + \epsilon_{\text{bend}} = \frac{w |\hat{\textbf{b}} \cdot \boldsymbol{\kappa}|}{2}, where :math:`\tau` is the torsion of the tape frame, :math:`\hat{\textbf{b}}` is the - frame binormal vector, and :math:`\boldsymbol{\kappa}` is the curvature vector of the - filamentary coil. + frame binormal vector, :math:`\boldsymbol{\kappa}` is the curvature vector of the + filamentary coil, and :math:`w` is the width of the tape. + This class is not intended to be used as an objective function inside + optimization. For that purpose you should instead use + :obj:`LPBinormalCurvatureStrainPenalty` or :obj:`LPTorsionalStrainPenalty`. + Those classes also compute gradients whereas this class does not. """ def __init__(self, framedcurve, width=1e-3): @@ -143,14 +161,14 @@ def __init__(self, framedcurve, width=1e-3): super().__init__(depends_on=[framedcurve]) def torsional_strain(self): - """ + r""" Returns the value of the torsional strain, :math:`\epsilon_{\text{tor}}`, along the quadpoints defining the filamentary coil. """ return self.torstrain_jax(self.framedcurve.frame_torsion(), self.width) def binormal_curvature_strain(self): - """ + r""" Returns the value of the torsional strain, :math:`\epsilon_{\text{bend}}`, along the quadpoints defining the filamentary coil. """ diff --git a/tests/geo/test_strainopt.py b/tests/geo/test_strainopt.py index 277aae347..87069508e 100644 --- a/tests/geo/test_strainopt.py +++ b/tests/geo/test_strainopt.py @@ -7,7 +7,7 @@ from scipy.optimize import minimize -class StrainOptTesting(unittest.TestCase): +class CoilStrainTesting(unittest.TestCase): def test_strain_opt(self): """