diff --git a/dev/index.html b/dev/index.html index fc7c196..8a17b53 100644 --- a/dev/index.html +++ b/dev/index.html @@ -6,4 +6,4 @@ & lcon \leq c(κ,y,u) \leq ucon,\\ & lvar \leq (κ,y,u) \leq uvar.\\ \end{aligned} -\right.\]

The main challenges in modeling such a problem are to be able to discretize the domain and generate corresponding discretizations of the objective and constraints, and their evaluate derivatives with respect to all variables. We use Gridap.jl to define the domain, meshes, function spaces, and finite-element families to approximate unknowns, and to model functionals and sets of PDEs in a weak form. PDENLPModels extends Gridap.jl's differentiation facilities to also obtain derivatives useful for optimization, i.e., first and second derivatives of the objective and constraint functions with respect to controls and finite-dimensional variables.

After discretization of the domain $\Omega$, the integral, and the derivatives, the resulting problem is a nonlinear optimization problem. PDENLPModels exports the GridapPDENLPModel type, an instance of an AbstractNLPModel, as defined in NLPModels.jl, which provides access to objective and constraint function values, to their first and second derivatives, and to any information that a solver might request from a model. The role of NLPModels.jl is to define an API that users and solvers can rely on. It is the role of other packages to implement facilities that create models compliant with the NLPModels API. We refer to juliasmoothoptimizers.github.io for tutorials on the NLPModel API.

As such, PDENLPModels offers an interface between generic PDE-constrained optimization problems and cutting-edge optimization solvers such as Artelys Knitro via NLPModelsKnitro.jl, Ipopt via NLPModelsIpopt.jl , DCISolver.jl, Percival.jl, and any solver accepting an AbstractNLPModel as input, see JuliaSmoothOptimizers.

Migot, T., Orban D., & Siqueira A. S. PDENLPModels.jl: A NLPModel API for optimization problems with PDE-constraints Journal of Open Source Software 7(80), 4736 (2022). 10.21105/joss.04736

Installation

] add PDENLPModels

The current version of PDENLPModels relies on Gridap v0.15.5.

Table of Contents

Examples

You can also check the tutorial Solve a PDE-constrained optimization problem on our site, juliasmoothoptimizers.github.io.

We refer to the folder test/problems for more examples of problems of different types: calculus of variations, optimal control problem, PDE-constrained problems, and mixed PDE-contrained problems with both function and vector unknowns. An alternative is to visit the repository PDEOptimizationProblems that contains a collection of test problems. Without objective function, the problem reduces to a classical PDE and we refer to Gridap tutorials for examples.

References

Gridap.jl Badia, S., Verdugo, F. (2020). Gridap: An extensible Finite Element toolbox in Julia. Journal of Open Source Software, 5(52), 2520.

NLPModels.jl D. Orban, A. S. Siqueira and contributors (2020). NLPModels.jl: Data Structures for Optimization Models

Bug reports and discussions

If you think you found a bug, feel free to open an issue. Focused suggestions and requests can also be opened as issues. Before opening a pull request, start an issue or a discussion on the topic, please.

If you want to ask a question not suited for a bug report, feel free to start a discussion here. This forum is for general discussion about this repository and the JuliaSmoothOptimizers, so questions about any of our packages are welcome.

+\right.\]

The main challenges in modeling such a problem are to be able to discretize the domain and generate corresponding discretizations of the objective and constraints, and their evaluate derivatives with respect to all variables. We use Gridap.jl to define the domain, meshes, function spaces, and finite-element families to approximate unknowns, and to model functionals and sets of PDEs in a weak form. PDENLPModels extends Gridap.jl's differentiation facilities to also obtain derivatives useful for optimization, i.e., first and second derivatives of the objective and constraint functions with respect to controls and finite-dimensional variables.

After discretization of the domain $\Omega$, the integral, and the derivatives, the resulting problem is a nonlinear optimization problem. PDENLPModels exports the GridapPDENLPModel type, an instance of an AbstractNLPModel, as defined in NLPModels.jl, which provides access to objective and constraint function values, to their first and second derivatives, and to any information that a solver might request from a model. The role of NLPModels.jl is to define an API that users and solvers can rely on. It is the role of other packages to implement facilities that create models compliant with the NLPModels API. We refer to juliasmoothoptimizers.github.io for tutorials on the NLPModel API.

As such, PDENLPModels offers an interface between generic PDE-constrained optimization problems and cutting-edge optimization solvers such as Artelys Knitro via NLPModelsKnitro.jl, Ipopt via NLPModelsIpopt.jl , DCISolver.jl, Percival.jl, and any solver accepting an AbstractNLPModel as input, see JuliaSmoothOptimizers.

Migot, T., Orban D., & Siqueira A. S. PDENLPModels.jl: A NLPModel API for optimization problems with PDE-constraints Journal of Open Source Software 7(80), 4736 (2022). 10.21105/joss.04736

Installation

] add PDENLPModels

The current version of PDENLPModels relies on Gridap v0.15.5.

Table of Contents

Examples

You can also check the tutorial Solve a PDE-constrained optimization problem on our site, juliasmoothoptimizers.github.io.

We refer to the folder test/problems for more examples of problems of different types: calculus of variations, optimal control problem, PDE-constrained problems, and mixed PDE-contrained problems with both function and vector unknowns. An alternative is to visit the repository PDEOptimizationProblems that contains a collection of test problems. Without objective function, the problem reduces to a classical PDE and we refer to Gridap tutorials for examples.

References

Gridap.jl Badia, S., Verdugo, F. (2020). Gridap: An extensible Finite Element toolbox in Julia. Journal of Open Source Software, 5(52), 2520.

NLPModels.jl D. Orban, A. S. Siqueira and contributors (2020). NLPModels.jl: Data Structures for Optimization Models

Bug reports and discussions

If you think you found a bug, feel free to open an issue. Focused suggestions and requests can also be opened as issues. Before opening a pull request, start an issue or a discussion on the topic, please.

If you want to ask a question not suited for a bug report, feel free to start a discussion here. This forum is for general discussion about this repository and the JuliaSmoothOptimizers, so questions about any of our packages are welcome.

diff --git a/dev/poisson-boltzman/index.html b/dev/poisson-boltzman/index.html index 311ebcc..edfc5c4 100644 --- a/dev/poisson-boltzman/index.html +++ b/dev/poisson-boltzman/index.html @@ -54,4 +54,4 @@ yh = FEFunction(nlp.pdemeta.Ypde, yfv) ufv = stats.solution[1+Gridap.FESpaces.num_free_dofs(nlp.pdemeta.Ypde):end] uh = FEFunction(nlp.pdemeta.Ycon, ufv) -writevtk(nlp.pdemeta.tnrj.trian,"results",cellfields=["uh"=>uh, "yh"=>yh])

Finally, the solution is obtained using any software reading VTK, e.g. Paraview.

Solution of P-B equationControl of P-B equation

References

[1] Michael J. Holst The Poisson-Boltzmann equation: Analysis and multilevel numerical solution. Applied Mathematics and CRPC, California Institute of Technology. (1994). Hols94d.pdf

+writevtk(nlp.pdemeta.tnrj.trian,"results",cellfields=["uh"=>uh, "yh"=>yh])

Finally, the solution is obtained using any software reading VTK, e.g. Paraview.

Solution of P-B equationControl of P-B equation

References

[1] Michael J. Holst The Poisson-Boltzmann equation: Analysis and multilevel numerical solution. Applied Mathematics and CRPC, California Institute of Technology. (1994). Hols94d.pdf

diff --git a/dev/reference/index.html b/dev/reference/index.html index 3acf46d..1f42d0f 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -1,7 +1,7 @@ Reference · PDENLPModels.jl
Edit on GitHub

Reference

​

Contents

​

​

Index

​

​

PDENLPModels.EnergyFETerm — Type

AbstractEnergyTerm modeling the objective function of the optimization problem

\[\begin{aligned} \int_{\Omega} f(y,u) d\Omega. -\end{aligned}\]

Constructor:

EnergyFETerm(:: Function, :: Triangulation, :: Measure)

See also: MixedEnergyFETerm, NoFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.GridapPDENLPModel — Type

GridapPDENLPModel returns an instance of an AbstractNLPModel using Gridap.jl for the discretization of the domain with finite-elements. Given a domain Ω ⊂ ℜᵈ Find a state function y: Ω -> Y, a control function u: Ω -> U and an algebraic vector κ ∈ ℜⁿ satisfying

\[\begin{aligned} +\end{aligned}\]

Constructor:

EnergyFETerm(:: Function, :: Triangulation, :: Measure)

See also: MixedEnergyFETerm, NoFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.GridapPDENLPModel — Type

GridapPDENLPModel returns an instance of an AbstractNLPModel using Gridap.jl for the discretization of the domain with finite-elements. Given a domain Ω ⊂ ℜᵈ Find a state function y: Ω -> Y, a control function u: Ω -> U and an algebraic vector κ ∈ ℜⁿ satisfying

\[\begin{aligned} \min_{κ,y,u} \ & ∫_Ω​ f(κ,y,u) dΩ​ \\ \mbox{ s.t. } & y \mbox{ solution of } PDE(κ,u)=0, \\ & lcon <= c(κ,y,u) <= ucon, \\ @@ -49,25 +49,25 @@ ncon = num_free_dofs(Ycon) xin = zeros(npde + ncon) - nlp = GridapPDENLPModel(xin, f, trian, Ypde, Ycon, Xpde, Xcon, res, name = "Control elastic membrane")

You can also check the tutorial Solve a PDE-constrained optimization problem on JSO's website, juliasmoothoptimizers.github.io.

We refer to the folder test/problems for more examples of problems of different types:

  • calculus of variations,
  • optimal control problem,
  • PDE-constrained problems,
  • mixed PDE-contrained problems with both function and algebraic unknowns.

An alternative is to visit the repository PDEOptimizationProblems that contains a collection of test problems.

Without objective function, the problem reduces to a classical PDE and we refer to Gridap tutorials for examples.

source
PDENLPModels.MixedEnergyFETerm — Type

AbstractEnergyTerm modeling the objective function of the optimization problem with functional and discrete unknowns

\[\begin{aligned} + nlp = GridapPDENLPModel(xin, f, trian, Ypde, Ycon, Xpde, Xcon, res, name = "Control elastic membrane")

You can also check the tutorial Solve a PDE-constrained optimization problem on JSO's website, juliasmoothoptimizers.github.io.

We refer to the folder test/problems for more examples of problems of different types:

  • calculus of variations,
  • optimal control problem,
  • PDE-constrained problems,
  • mixed PDE-contrained problems with both function and algebraic unknowns.

An alternative is to visit the repository PDEOptimizationProblems that contains a collection of test problems.

Without objective function, the problem reduces to a classical PDE and we refer to Gridap tutorials for examples.

source
PDENLPModels.MixedEnergyFETerm — Type

AbstractEnergyTerm modeling the objective function of the optimization problem with functional and discrete unknowns

\[\begin{aligned} \int_{\Omega} f(y,u,\kappa) d\Omega. -\end{aligned}\]

Constructor:

MixedEnergyFETerm(:: Function, :: Triangulation, :: Int)

See also: EnergyFETerm, NoFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.NoFETerm — Type

AbstractEnergyTerm modeling the objective function when there are no integral objective

\[\begin{aligned} +\end{aligned}\]

Constructor:

MixedEnergyFETerm(:: Function, :: Triangulation, :: Int)

See also: EnergyFETerm, NoFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.NoFETerm — Type

AbstractEnergyTerm modeling the objective function when there are no integral objective

\[\begin{aligned} f(\kappa). -\end{aligned}\]

Constructors: NoFETerm() NoFETerm(:: Function)

See also: MixedEnergyFETerm, EnergyFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.PDENLPMeta — Type
PDENLPMeta

A composite type that represents the main features of the PDE-constrained optimization problem

PDENLPMeta contains the following attributes:

  • tnrj : structure representing the objective term
  • Ypde : TrialFESpace for the solution of the PDE
  • Ycon : TrialFESpace for the parameter
  • Xpde : TestFESpace for the solution of the PDE
  • Xcon : TestFESpace for the parameter
  • Y : concatenated TrialFESpace
  • X : concatenated TestFESpace
  • op : operator representing the PDE-constraint (nothing if no constraints)
  • nvar_pde :number of dofs in the solution functions
  • nvar_con : number of dofs in the control functions
  • nparam : number of real unknowns
  • nnzh_obj : number of nonzeros elements in the objective hessian
  • Hrows : store the structure for the hessian of the lagrangian
  • Hcols : store the structure for the hessian of the lagrangian
  • Jrows : store the structure for the hessian of the jacobian
  • Jcols : store the structure for the hessian of the jacobian
source
PDENLPModels._compute_gradient! — Function

Return the gradient of the objective function and set it in place.

_compute_gradient!(:: AbstractVector, :: EnergyFETerm, :: AbstractVector, :: FEFunctionType, :: FESpace, :: FESpace)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_hess_k_coo

source
PDENLPModels._compute_gradient_k! — Function

Return the derivative of the objective function w.r.t. κ.

_compute_gradient_k!(::AbstractVector, :: AbstractEnergyTerm, :: FEFunctionType, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_hess_k_coo

source
PDENLPModels._compute_hess_k_vals! — Function

Return the values of the hessian w.r.t. κ of the objective function.

_compute_hess_k_vals!(:: AbstractVector, :: AbstractNLPModel, :: AbstractEnergyTerm, :: AbstractVector, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_gradient_k

source
PDENLPModels._compute_hess_structure — Method
_compute_hess_structure(::AbstractEnergyTerm, Y, X, x0, nparam)
-_compute_hess_structure::AbstractEnergyTerm, op, Y, Ypde, Ycon, X, Xpde, x0, nparam)

Return a triplet with the structure (rows, cols) and the number of non-zeros elements in the hessian w.r.t. y and u of the objective function.

The rows and cols returned by _compute_hess_structure_obj are already shifter by nparam.

source
PDENLPModels._functions_to_vectors! — Method
_functions_to_vectors!(nini :: Int, nfields :: Int, trian :: Triangulation, lfunc :: Function, ufunc :: Function, cell_xm :: Gridap.Arrays.AppliedArray, Y :: FESpace, lvar :: AbstractVector, uvar :: AbstractVector)

Iterate for k = 1 to nfields and switch lfunc[k] and ufunc[k] to vectors, allocated in lvar and uvar in place starting from nini + 1. It returns nini + the number of allocations.

source
PDENLPModels._obj_integral — Function

Return the integral of the objective function

_obj_integral(:: AbstractEnergyTerm, :: FEFunctionType, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_cell_integral, _compute_gradient_k, _compute_hess_k_coo

source
PDENLPModels._split_FEFunction — Method
_split_FEFunction(:: AbstractVector,  :: FESpace, :: Union{FESpace, Nothing})

Split the vector x into two FEFunction corresponding to the solution y and the control u. Returns nothing for the control u if Ycon == nothing.

Do not verify the compatible length.

source
PDENLPModels._split_vector — Method
_split_vector(:: AbstractVector,  :: FESpace, :: Union{FESpace, Nothing})

Split the vector x into three vectors: y, u, k. Returns nothing for the control u if Ycon == nothing.

Do not verify the compatible length.

source
PDENLPModels.bounds_functions_to_vectors — Method
bounds_functions_to_vectors(Y :: MultiFieldFESpace, Ycon :: Union{FESpace, Nothing},  Ypde :: FESpace, trian :: Triangulation, lyfunc :: Union{Function, AbstractVector}, uyfunc :: Union{Function, AbstractVector}, lufunc :: Union{Function, AbstractVector}, uufunc :: Union{Function, AbstractVector})

Return the bounds lvar and uvar.

source
+\end{aligned}\]

Constructors: NoFETerm() NoFETerm(:: Function)

See also: MixedEnergyFETerm, EnergyFETerm, _obj_cell_integral, _obj_integral, _compute_gradient_k!

source
PDENLPModels.PDENLPMeta — Type
PDENLPMeta

A composite type that represents the main features of the PDE-constrained optimization problem

PDENLPMeta contains the following attributes:

  • tnrj : structure representing the objective term
  • Ypde : TrialFESpace for the solution of the PDE
  • Ycon : TrialFESpace for the parameter
  • Xpde : TestFESpace for the solution of the PDE
  • Xcon : TestFESpace for the parameter
  • Y : concatenated TrialFESpace
  • X : concatenated TestFESpace
  • op : operator representing the PDE-constraint (nothing if no constraints)
  • nvar_pde :number of dofs in the solution functions
  • nvar_con : number of dofs in the control functions
  • nparam : number of real unknowns
  • nnzh_obj : number of nonzeros elements in the objective hessian
  • Hrows : store the structure for the hessian of the lagrangian
  • Hcols : store the structure for the hessian of the lagrangian
  • Jrows : store the structure for the hessian of the jacobian
  • Jcols : store the structure for the hessian of the jacobian
source
PDENLPModels.PDEWorkspace — Type
PDEWorkspace

Pre-allocated memory for GridapPDENLPModel.

source
PDENLPModels._compute_gradient! — Function

Return the gradient of the objective function and set it in place.

_compute_gradient!(:: AbstractVector, :: EnergyFETerm, :: AbstractVector, :: FEFunctionType, :: FESpace, :: FESpace)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_hess_k_coo

source
PDENLPModels._compute_gradient_k! — Function

Return the derivative of the objective function w.r.t. κ.

_compute_gradient_k!(::AbstractVector, :: AbstractEnergyTerm, :: FEFunctionType, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_hess_k_coo

source
PDENLPModels._compute_hess_k_vals! — Function

Return the values of the hessian w.r.t. κ of the objective function.

_compute_hess_k_vals!(:: AbstractVector, :: AbstractNLPModel, :: AbstractEnergyTerm, :: AbstractVector, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_integral, _obj_cell_integral, _compute_gradient_k

source
PDENLPModels._compute_hess_structure — Method
_compute_hess_structure(::AbstractEnergyTerm, Y, X, x0, nparam)
+_compute_hess_structure::AbstractEnergyTerm, op, Y, Ypde, Ycon, X, Xpde, x0, nparam)

Return a triplet with the structure (rows, cols) and the number of non-zeros elements in the hessian w.r.t. y and u of the objective function.

The rows and cols returned by _compute_hess_structure_obj are already shifter by nparam.

source
PDENLPModels._functions_to_vectors! — Method
_functions_to_vectors!(nini :: Int, nfields :: Int, trian :: Triangulation, lfunc :: Function, ufunc :: Function, cell_xm :: Gridap.Arrays.AppliedArray, Y :: FESpace, lvar :: AbstractVector, uvar :: AbstractVector)

Iterate for k = 1 to nfields and switch lfunc[k] and ufunc[k] to vectors, allocated in lvar and uvar in place starting from nini + 1. It returns nini + the number of allocations.

source
PDENLPModels._obj_integral — Function

Return the integral of the objective function

_obj_integral(:: AbstractEnergyTerm, :: FEFunctionType, :: AbstractVector)

See also: MixedEnergyFETerm, EnergyFETerm, NoFETerm, _obj_cell_integral, _compute_gradient_k, _compute_hess_k_coo

source
PDENLPModels._split_FEFunction — Method
_split_FEFunction(:: AbstractVector,  :: FESpace, :: Union{FESpace, Nothing})

Split the vector x into two FEFunction corresponding to the solution y and the control u. Returns nothing for the control u if Ycon == nothing.

Do not verify the compatible length.

source
PDENLPModels._split_vector — Method
_split_vector(:: AbstractVector,  :: FESpace, :: Union{FESpace, Nothing})

Split the vector x into three vectors: y, u, k. Returns nothing for the control u if Ycon == nothing.

Do not verify the compatible length.

source
PDENLPModels._split_vectors — Method
_split_vectors(x, Ypde, Ycon)

Take a vector x and returns a splitting in terms of y, u and θ.

source
PDENLPModels.bounds_functions_to_vectors — Method
bounds_functions_to_vectors(Y :: MultiFieldFESpace, Ycon :: Union{FESpace, Nothing},  Ypde :: FESpace, trian :: Triangulation, lyfunc :: Union{Function, AbstractVector}, uyfunc :: Union{Function, AbstractVector}, lufunc :: Union{Function, AbstractVector}, uufunc :: Union{Function, AbstractVector})

Return the bounds lvar and uvar.

source
PDENLPModels.get_Hcols — Method
get_Hcols(nlp)
+get_Hcols(meta)

Return the value Hcols from meta or nlp.meta.

source
PDENLPModels.get_Hrows — Method
get_Hrows(nlp)
+get_Hrows(meta)

Return the value Hrows from meta or nlp.meta.

source
PDENLPModels.get_Jcols — Method
get_Jcols(nlp)
+get_Jcols(meta)

Return the value Jcols from meta or nlp.meta.

source
PDENLPModels.get_Jrows — Method
get_Jrows(nlp)
+get_Jrows(meta)

Return the value Jrows from meta or nlp.meta.

source
PDENLPModels.get_X — Method
get_X(nlp)
+get_X(meta)

Return the value X from meta or nlp.meta.

source
PDENLPModels.get_Xcon — Method
get_Xcon(nlp)
+get_Xcon(meta)

Return the value Xcon from meta or nlp.meta.

source
PDENLPModels.get_Xpde — Method
get_Xpde(nlp)
+get_Xpde(meta)

Return the value Xpde from meta or nlp.meta.

source
PDENLPModels.get_Y — Method
get_Y(nlp)
+get_Y(meta)

Return the value Y from meta or nlp.meta.

source
PDENLPModels.get_Ycon — Method
get_Ycon(nlp)
+get_Ycon(meta)

Return the value Ycon from meta or nlp.meta.

source
PDENLPModels.get_Ypde — Method
get_Ypde(nlp)
+get_Ypde(meta)

Return the value Ypde from meta or nlp.meta.

source
PDENLPModels.get_nnzh_obj — Method
get_nnzh_obj(nlp)
+get_nnzh_obj(meta)

Return the value nnzh_obj from meta or nlp.meta.

source
PDENLPModels.get_nparam — Method
get_nparam(nlp)
+get_nparam(meta)

Return the value nparam from meta or nlp.meta.

source
PDENLPModels.get_nvar_con — Method
get_nvar_con(nlp)
+get_nvar_con(meta)

Return the value nvar_con from meta or nlp.meta.

source
PDENLPModels.get_nvar_pde — Method
get_nvar_pde(nlp)
+get_nvar_pde(meta)

Return the value nvar_pde from meta or nlp.meta.

source
PDENLPModels.get_op — Method
get_op(nlp)
+get_op(meta)

Return the value op from meta or nlp.meta.

source
PDENLPModels.get_tnrj — Method
get_tnrj(nlp)
+get_tnrj(meta)

Return the value tnrj from meta or nlp.meta.

source
PDENLPModels.split_vectors — Method
split_vectors(::GridapPDENLPModel, x)

Take a vector x and returns a splitting in terms of y, u and θ.

source
diff --git a/dev/search/index.html b/dev/search/index.html index e4fa8cb..4ea0722 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -Search · PDENLPModels.jl

Loading search...

    +Search · PDENLPModels.jl

    Loading search...

      diff --git a/dev/tore/index.html b/dev/tore/index.html index 62fb43f..216fa24 100644 --- a/dev/tore/index.html +++ b/dev/tore/index.html @@ -96,4 +96,4 @@ Ys = (c .+ a * cos.(αs)) * sin.(βs)' Zs = (a * sin.(αs)) * ones(M)' plot3d(Xs, Ys, Zs, st=:surface, grid=false, c=:grays, axis=false, colorbar=false) -plot3d!(xs, ys, zs, linewidth=4, color=:red, title=@sprintf("Geodesic on a Torus (length=%4.f)", L), legend=false)

      Geodesic over the tore

      References

      [1] Weisstein, Eric W. Brachistochrone Problem. From MathWorld–A Wolfram Web Resource. mathworld.wolfram.com/BrachistochroneProblem.html

      [2] Mark L. Irons The Curvature and Geodesics of the Torus (2005). torus.geodesics.pdf

      +plot3d!(xs, ys, zs, linewidth=4, color=:red, title=@sprintf("Geodesic on a Torus (length=%4.f)", L), legend=false)

      Geodesic over the tore

      References

      [1] Weisstein, Eric W. Brachistochrone Problem. From MathWorld–A Wolfram Web Resource. mathworld.wolfram.com/BrachistochroneProblem.html

      [2] Mark L. Irons The Curvature and Geodesics of the Torus (2005). torus.geodesics.pdf