Giter Club home page Giter Club logo

pacopy's Issues

Add very simple examples for beginners

Hi,

thank you for this great module! I think the included examples are fairly complicated, wouldn't it be nice to have some very simple examples? In my opinion a simple sin(u)-lmbda would be enough to show the "curve following process" of the given equilibrium equation(s).

My suggestion would be something like:

#from autograd import numpy
#import autograd
import numpy
import scipy
import matplotlib.pyplot as plt

import pacopy

class SimpleProblem1D:
    def __init__(self):
        self.n = 1

    def inner(self, a, b):
        """The inner product of the problem. Can be numpy.dot(a, b), but factoring in
        the mesh width stays true to the PDE.
        """
        return numpy.dot(a, b)

    def norm2_r(self, a):
        """The norm in the range space; used to determine if a solution has been found."""
        return numpy.dot(a, a)

    def f(self, u, lmbda):
        """The evaluation of the function to be solved"""
        out = numpy.sin(u) - lmbda
        return out

    def df_dlmbda(self, u, lmbda):
        """The function's derivative with respect to the parameter. Used in Euler-Newton
        continuation.
        """
        out = -numpy.ones_like(u)
        return out

    def jacobian_solver(self, u, lmbda, rhs):
        """A solver for the Jacobian problem."""
        #df_du = autograd.jacobian(self.f)(u, lmbda)
        df_du = numpy.array([numpy.cos(u)])
        return numpy.linalg.solve(df_du, rhs)
    
problem = SimpleProblem1D()
# Initial guess
u0 = numpy.zeros(problem.n)
# Initial parameter value
lmbda0 = 0.0

# Init lists
lmbda_list = []
values_list = []

def callback(k, lmbda, sol):
    # Use the callback for plotting, writing data to files etc.
    lmbda_list.append(lmbda)
    values_list.append(sol)

pacopy.euler_newton(problem, u0, lmbda0, callback, max_steps=20, newton_tol=1.0e-10)

# plot solution
fig = plt.figure()
data = numpy.hstack((numpy.array(values_list),
                     numpy.array(lmbda_list).reshape(-1,1))).T
plt.plot(*data,'.-')
plt.xlabel('$u_1$')
plt.ylabel('$\lambda$')

Best regards,
Andreas

no code in main branch

Hi,
I'm wondering why there is no code in the main branch?
If I'd like to use pacopy shall I consider using branch-switching instead

Best

version bump for pacopy.natural.milestones

pacopy.natural.milestones #5 is working nicely but I just realized that it's not yet available with pip install pacopy which is still at v. 0.1.0.

I'd submit a PR but I'm not sure whether it would mean 0.1.1 or 0.2.0.

milestones

Would it be useful to add an optional argument to natural and euler_newton to force solution at certain values of λ ?

It's easy enough now to reproduce, as in test/test_bratu1d.py, most of the left plot in figure 1.1 of ‘Deflation techniques…’ (Farrell, Birkisson, & Funke), but not the blue or red cross at λ = 2 or the right plot with the corresponding solutions u (x; λ = 2).

scipy.integrate.solve_ivp has a couple of related features:

  • t_eval: ‘Times at which to store the computed solution’
  • events: ‘defined by a continuous function of time and state that becomes zero value in case of an event’

I envisage something like replacing

https://github.com/nschloe/pacopy/blob/1973500e422132f035fd9a419340c5100d6bb904/pacopy/natural.py#L79

with

milestone = next(milestones)
...
while ...:
    lmbda_next = lmbda + lmbda_stepsize
    if milestone < lmbda_next:
        lmbda = milestone
        milestone = next(milestones)
    else:
        lmbda = lmbda_next

plus handing:

  • lmbda_stepsize might be negative
  • milestones might be empty

The motivation is a steady Navier–Stokes solver for which the solutions are desired at certain Reynolds numbers.

This version of pacopy is outdated. Please upgrade.

I get this error even though I download with latest pip not sure what is going on.
Successfully installed markdown-it-py-3.0.0 mdurl-0.1.2 pacopy-0.2.7

pacopy-0.2.7 is the same one available on the main page of this repo at the moment

default maximum parameter stepsize

The default maximum parameter stepsize is finite

https://github.com/nschloe/pacopy/blob/b5ae74c2c17349a04165291e2e74fbc865da660b/pacopy/natural.py#L12

and small, depending on the system. A couple of times I have forgotten this and removed the explicit setting, hoping for a faster run through the milestones #6 (Δλ ~ 10², in the example given there, kinnala/scikit-fem#145), only to have it tiptoe at Δλ = 1.

What about defaulting to no limit as for

https://github.com/nschloe/pacopy/blob/b5ae74c2c17349a04165291e2e74fbc865da660b/pacopy/natural.py#L16

?

(Otherwise, pacopy.natural is proving very useful for stationary Navier–Stokes work.)

readme dump

This thread is for storing videos to be used in the readme.

Complex-Ginsburg-Landau

I'm getting this broadcasting error for the shapes and having trouble sourcing the error.... any idea? Thanks!

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (3,19602,3)->(3,19602,3) (3,19602,2)->(3,19602,2)"

This is the error that is returned when code is run:

  # test_f_i_psi()
   # test_df_dlmbda()
--->      test_ginzburg_landau(max_steps=100, n=100)
   # plot_data()'''

<ipython-input-1-cd6944349cb4> in test_ginzburg_landau(max_steps, n)'''
   302         #     newton_tol=1.0e-10,
   303         # )
--> 304         pacopy.euler_newton(
   305             problem,
   306             u0,

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/pacopy/euler_newton.py in euler_newton(problem, u0, lmbda0, callback, max_steps, verbose, newton_tol, max_newton_steps, predictor_variant, corrector_variant, stepsize0, stepsize_max, stepsize_aggressiveness, cos_alpha_min, theta0, adaptive_theta, converge_onto_zero_eigenvalue)'''
   140     k = 0
   141     try:
--> 142         u, _ = newton(
   143             lambda u: problem.f(u, lmbda),
   144             lambda u, rhs: problem.jacobian_solver(u, lmbda, rhs),

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/pacopy/newton.py in newton(f, jacobian_solver, norm2, u0, tol, max_iter, verbose)'''
     9     u = u0
    10 
---> 11     fu = f(u)
    12     nrm = math.sqrt(norm2(fu))
    13     if verbose:

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/pacopy/euler_newton.py in <lambda>(u)'''
   141     try:
   142         u, _ = newton(
--> 143             lambda u: problem.f(u, lmbda),
   144             lambda u, rhs: problem.jacobian_solver(u, lmbda, rhs),
   145             problem.norm2_r,

<ipython-input-1-cd6944349cb4> in f(self, psi, mu)'''
    99     def f(self, psi, mu):
--> 100         keo = pyfvm.get_fvm_matrix(self.mesh, edge_kernels=[Energy(mu)])
   101         cv = self.mesh.control_volumes
   102         out = (keo * psi) / cv + (self.V + self.g * numpy.abs(psi) ** 2) * psi

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/pyfvm/fvm_matrix.py in get_fvm_matrix(mesh, edge_kernels, vertex_kernels, face_kernels, dirichlets)'''
    11     dirichlets = [] if dirichlets is None else dirichlets
    12 
---> 13     V, I, J = _get_VIJ(mesh, edge_kernels, vertex_kernels, face_kernels)
    14 
    15     # One unknown per vertex

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/pyfvm/fvm_matrix.py in _get_VIJ(mesh, edge_kernels, vertex_kernels, face_kernels)'''
    44             cell_mask = mesh.get_cell_mask(subdomain)
    45 
---> 46             v_matrix = edge_kernel.eval(mesh, cell_mask)
    47 
    48             V.append(v_matrix[0, 0].flatten())

<ipython-input-1-cd6944349cb4> in eval(self, mesh, cell_mask)'''
    35         # The dot product <magnetic_potential, edge>, executed for many
    36         # points at once; cf. <http://stackoverflow.com/a/26168677/353337>.
---> 37         beta = numpy.einsum("...k,...k->...", magnetic_potential, edge)
    38 
    39         return numpy.array( <__array_function__ internals> in einsum(*args, **kwargs)

~/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/numpy/core/einsumfunc.py in einsum(out, optimize, *operands, **kwargs)'''
  1348         if specified_out:
  1349             kwargs['out'] = out
-> 1350         return c_einsum(*operands, **kwargs)
  1351 
  1352     # Check the kwargs to avoid a more cryptic error later, without having to

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (3,19602,3)->(3,19602,3) (3,19602,2)->(3,19602,2) 

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.