diff --git a/doc/whats-new.md b/doc/whats-new.md index b2e376a4..593b91a2 100644 --- a/doc/whats-new.md +++ b/doc/whats-new.md @@ -25,6 +25,9 @@ Release history - Wall coordinates are written to output grid as `closed_wall_R` and `closed_wall_Z` (#176) By [Ben Dudson](https://github.com/bendudson) +- Extend divertor legs with `leg_extend` options. These specify how far each + leg should extend beyond the wall intersection (#195). + By [Ben Dudson](https://github.com/bendudson) 0.5.2 (13th March 2023) ------------------------- diff --git a/hypnotoad/cases/tokamak.py b/hypnotoad/cases/tokamak.py index 4a472c31..86c52bf0 100644 --- a/hypnotoad/cases/tokamak.py +++ b/hypnotoad/cases/tokamak.py @@ -78,6 +78,36 @@ class TokamakEquilibrium(Equilibrium): value_type=float, check_all=is_positive, ), + leg_extend=WithMeta( + 0.0, + doc="Extend divertor legs by this length [m]", + value_type=float, + check_all=is_non_negative, + ), + leg_extend_lower_inner=WithMeta( + "leg_extend", + doc="Extend lower inner divertor leg", + value_type=float, + check_all=is_non_negative, + ), + leg_extend_lower_outer=WithMeta( + "leg_extend", + doc="Extend lower outer divertor leg", + value_type=float, + check_all=is_non_negative, + ), + leg_extend_upper_inner=WithMeta( + "leg_extend", + doc="Extend upper inner divertor leg", + value_type=float, + check_all=is_non_negative, + ), + leg_extend_upper_outer=WithMeta( + "leg_extend", + doc="Extend upper outer divertor leg", + value_type=float, + check_all=is_non_negative, + ), nx_core=WithMeta( 5, doc="Number of radial points in the core", @@ -616,6 +646,20 @@ def findLegs(self, xpoint, radius=0.01, step=0.01): # The sign is used to tell which way to integrate sign = np.sign((leg[0] - xpoint.R) * Br + (leg[1] - xpoint.Z) * Bz) + # Get the length that the leg should be extended by + if leg[1] > xpoint.Z: + # Upper + if leg[0] > xpoint.R: + leg_extend = self.user_options.leg_extend_upper_outer + else: + leg_extend = self.user_options.leg_extend_upper_inner + else: + # Lower + if leg[0] > xpoint.R: + leg_extend = self.user_options.leg_extend_lower_outer + else: + leg_extend = self.user_options.leg_extend_lower_inner + # Integrate in this direction until the wall is intersected # This is affected by sign, which determines which way to integrate def dpos_dl(distance, pos): @@ -644,6 +688,20 @@ def dpos_dl(distance, pos): intersect = self.wallIntersection(Point2D(*pos), Point2D(*newpos)) if intersect is not None: line.append(intersect) # Put the intersection in the line + if leg_extend > 0.0: + nsteps = int(leg_extend / step + 0.5) + extend_step = leg_extend / nsteps + pos = (intersect.R, intersect.Z) + for i in range(nsteps): + solve_result = solve_ivp( + dpos_dl, + (0.0, extend_step), + pos, + rtol=0.0, + atol=self.user_options.leg_trace_atol, + ) + pos = (solve_result.y[0][1], solve_result.y[1][1]) + line.append(Point2D(*pos)) break pos = newpos line.append(Point2D(*pos))