Point evaluations #612

miklos1 · 2015-10-16T16:12:07Z

What is still missing:

Immersed manifolds: While much of the infrastructure already supports evaluating points on immersed manifolds, ~~finding the point on the manifold may or may not succeed.~~ this feature is disabled for now.

Limited support:

Distributed mesh support: Works, but parallel efficiency (to be honest: any efficiency) is not a concern at this stage. Point evaluation is a collective operation, and the same set of points are expected to be supplied on each process, and each process will have the values for each point.

Usage:

f.at(points) with all kinds of fancy syntax, e.g. f.at(1.2, 4.5) (one point), f.at([0.5, -2.1], [-2.6, 3.4]) (multiple points) or f.at(np.array([[0.5, -2.1], [-2.6, 3.4]])) (array of points). An exception is thrown if any point is outside the domain, unless f.at(points, dont_raise=True) is used, in which case the the values for points outside the domain will be filled with Nones.
UFL allows the evaluation of any expression at a single point with the function call syntax. For example, (f*sin(f))([0.3, 0.8]).

This depends on OP2/PyOP2#472 and on https://bitbucket.org/mapdes/ffc/pull-requests/86.

TODO: - Building C++

dorugeber · 2015-10-16T16:17:37Z

.travis.yml

@@ -49,7 +50,7 @@ before_install:
 install:
  - mkdir tmp
  - cd tmp
-  - ../scripts/firedrake-install --disable_ssh --global --minimal_petsc
+  - ../scripts/firedrake-install --disable_ssh --global --minimal_petsc --package_branch PyOP2 refactor-wrapper --package_branch ffc compile-element


STILL NEED TO DROP BEFORE MERGE

miklos1 · 2015-10-20T10:34:09Z

@dham: This will add a new dependency (libspatialindex). Will firedrake-update try to install the new dependency, or this feature does not exist yet? Or shall we just e-mail the Firedrake list so that the users will know to that they need to install the new dependency?

dham · 2015-10-20T10:40:25Z

You'll need to explicitly add the dependency to the apt and brew
dependencies in firedrake-install. This will probably entail users
rebuilding firedrake-update using the new instructions.

Cheers,

David

On Tue, 20 Oct 2015 at 11:34 Miklós Homolya [email protected]
wrote:

@dham https://github.com/dham: This will add a new dependency
(libspatialindex). Will firedrake-update try to install the new
dependency, or this feature does not exist yet? Or shall we just e-mail the
Firedrake list so that the users will know to that they need to install the
new dependency?

—
Reply to this email directly or view it on GitHub
#612 (comment)
.

miklos1 · 2015-10-20T10:46:05Z

@dham: So I did that, and that's why I'm not worried about fresh installs. My question is what happens when we merge this and someone just tries to naively update with firedrake-update.

dham · 2015-10-20T11:00:44Z

Potentially they get a build fail. They'd have to update firedrake update.

This would be more robust if firedrake-update updated itself first.

On Tue, 20 Oct 2015 at 11:46 Miklós Homolya [email protected]
wrote:

@dham https://github.com/dham: So I did that, and that's why I'm not
worried about fresh installs. My question is what happens when we merge
this and someone just tries to naively update with firedrake-update.

—
Reply to this email directly or view it on GitHub
#612 (comment)
.

miklos1 · 2015-10-21T16:40:19Z

So if one wants to evaluate multiple points on a mixed function, how shall the result look like?
a) an array of mixed values

array([[3, array([1, 2])],
       [4, array([5, 6])]], dtype=object)

b) a tuple of arrays

(array([3, 4]), array([[1, 2], [5, 6]]))

miklos1 · 2015-10-27T16:48:30Z

So I have prepared this branch for merge, and I can merge as soon as the dependent branches are merged / approved.

wence- · 2015-10-28T11:25:55Z

docs/source/point-evaluation.rst

+
+For a single point, the result is ``numpy`` array, or a tuple of
+``numpy`` arrays in case *mixed* functions.  When evaluating multiple
+points, the result is a list of values for each point.


I would have expected to get a numpy array for multiple points, rather than a list?

Also: "The result is a numpy array". For mixed functions, I presume that you get one numpy array per subspace?

Single point, non-mixed: numpy array
Single point, mixed: tuple of numpy arrays
Multiple points, non-mixed: list of numpy arrays
Multiple points, mixed: list of tuples of numpy arrays

We could have a "2D" numpy array for multiple points, non-mixed, but the 2D array idea doesn't quite compose with mixed. Unless maybe, if in the multiple points, mixed case we return a tuple of 2D numpy arrays, but David didn't like that.

Oh, I had thought that in the mixed case, I'd get the equivalent of:

tuple(f.at(pts) for f in mixed_f)

But instead I get:

[tuple(f.at(pt) for f in mixed) for pt in pts]
?

Yes, currently the second version happens. But I'm willing to do the first if @dham agrees.

miklos1 · 2015-10-28T15:05:41Z

I think I have addressed all issues so far.

* point-evaluation-serial: review fixes to test cases review fixes to code and documentation add section about point evaluation in the manual a parallel test make parallel evaluation work (inefficiently) check input for parallel consistency fix parallel evaluation bug drop users ability to specify fill value mixed function test cases and fixes fix evaluation of mixed functions disable broken immersed support deduplicate "firedrake" PointOutOfDomainError -> PointNotInDomainError support TensorFunctionSpace dissolve cfunction.py revise and test API refactor mesh.spatial_index spatialindex.pyx docstrings split spatialindex.pyx refactor locate.cpp n_cells -> n_cols make cells affine for the function space test cases [FFC] pass cdim to compile_element test case: function spaces extruded test case: function spaces non-extruded test case: point in cell query extruded test case: point in cell query non-extruded add libspatialindex dependency use fancy C++ kernel fall back on brute force for 1-dimension initial support for extruded spatial index for extruded meshes [FFC] get wrapper from PyOP2 convenient API for function evaluation use libspatialindex to speed up lookups merge header files starting point evaluations

* Change requirements on loopy to point at the loopy sprint branch. * Replace lp.register_callable_kernel with lp.merge knl = lp.register_callable_kernel(knl, callee) is now knl = lp.merge([knl, callee]) . Reason given for loopy breaking back-compat: earlier they had a notion of a special "root_kernel" which only permitted one level of nesting loopy kernel calls and that was too restrictive. * Loopy callables: adapt to new API for registering. * Loopy Callables: register callables on a kernel not functions. * name is static property on SolveCallable and INVCallable * Loopy callables: register the callables rather than functions and remove the lookup functions. * Some progress * Loopy callables: Adapt API of with_types. * INVCallable has name "inverse" in loopy rep * Always pass a target when generating loopy kernels, otherwise we cannot merge them with wrapper kernels. * flake8 * Maybe generate call to inner kernel correctly Still doesn't work later. * Revert "Maybe generate call to inner kernel correctly" This reverts commit 009ca19f3fb133d53affeab53cd016c9b55620a9. * Revert "Revert "Maybe generate call to inner kernel correctly"" This reverts commit c85b60e6382b4762076b2744227fa5a6d903b830. * Remove dead code * rep2loopy: get name of callee kernel for matching the args. * Add missing code target * Maybe handle loopy argument access correctly dat.zero still doesn't work * Promote dats with shape (n,) to shape (n, 1) This enables sweeping indices to pick up a direct loop. * Set lang_version on all make_function calls * Fix disjoint arguments to C-string kernels * Resolve review comments * Rename function_args -> kernel_parameters * WIP Match args for all kernels in the callables table besides for LACallables? Why? * Only match args for callables that are CallableKernels. * Pin decorator Decorator version 5 changes the API in ways that break the argument type checking in PyOP2. We're therefore pinning on the version that works. * catch up with numpy type changes * Set within_inames_is_final on instructions * Squash numpy warning * Set strides on GlobalArgs * Caching: produce key for loop Programs. * Check if kernel is a loopy Program. * loopy.program.Program -> loopy.Program * Adapt to loopy interface changes. * Maintain loopy automatic reshaping (was in loopy) Everything in loopycompat.py file was in loopy/transform/callable.py but is removed in inducer/loopy#327. This commit maintains compatibility but should be phased out. * Fix loopy.program -> loopy.translation_unit * dont lint automatic reshaping (was in loopy) * Lint to pass flake8 * Fix loopy automatic reshaping (was in loopy) Changes were copied from an out of date branch by accident. * Fix type case error (recursively) This mirrors inducer/loopy#326 * ? Maybe fix the last dimension error mismatch. Not sure this will break things in other places. * Codegen: Don't merge indices of extent 1, so that the instructions are schedulable. This will fix e.g. the failure of tests/extrusion/test_steady_advection_2D_extr.py * Testing: Update loopy branch * Squashed commit of the following: commit 5a89bfc Author: Connor Ward <[email protected]> Date: Mon May 10 10:44:33 2021 +0100 Use Firedrake PETSc and petsc4py commit 948cb7d Author: David Ham <[email protected]> Date: Wed Apr 21 13:30:50 2021 +0100 Pin decorator Decorator version 5 changes the API in ways that break the argument type checking in PyOP2. We're therefore pinning on the version that works. commit 121c2a6 Author: David Ham <[email protected]> Date: Wed Apr 21 13:45:58 2021 +0100 catch up with numpy type changes * Revert "Squashed commit of the following:" This reverts commit bd1916142f883d6a4f0371ca6de2c344f0c7e78a. * Squashed commit of the following: commit 9a90c9f Merge: f1daf21 c8d7ff6 Author: Lawrence Mitchell <[email protected]> Date: Tue May 11 12:07:33 2021 +0100 Merge pull request #618 from OP2/wence/fix/actions-doc-build actions: Build docs at end of successful test run commit c8d7ff6 Author: Lawrence Mitchell <[email protected]> Date: Tue May 11 11:41:20 2021 +0100 actions: Build docs at end of successful test run commit f1daf21 Merge: d0cc348 5a89bfc Author: Lawrence Mitchell <[email protected]> Date: Mon May 10 20:32:24 2021 +0100 Merge pull request #617 from connorjward/test-ci Make PETSc and petsc4py for CI use Firedrake not pip commit 5a89bfc Author: Connor Ward <[email protected]> Date: Mon May 10 10:44:33 2021 +0100 Use Firedrake PETSc and petsc4py commit 948cb7d Author: David Ham <[email protected]> Date: Wed Apr 21 13:30:50 2021 +0100 Pin decorator Decorator version 5 changes the API in ways that break the argument type checking in PyOP2. We're therefore pinning on the version that works. commit 121c2a6 Author: David Ham <[email protected]> Date: Wed Apr 21 13:45:58 2021 +0100 catch up with numpy type changes * Lint. * Jenkins * Make import global. * Drop package branches. Co-authored-by: Sophia Vorderwuelbecke <[email protected]> Co-authored-by: Lawrence Mitchell <[email protected]> Co-authored-by: Connor Ward <[email protected]> Co-authored-by: David Ham <[email protected]> Co-authored-by: Kaushik Kulkarni <[email protected]>

miklos1 added 22 commits October 15, 2015 16:53

starting point evaluations

1b8a56d

merge header files

052f3e0

use libspatialindex to speed up lookups

140b6be

TODO: - Building C++

convenient API for function evaluation

a280237

[FFC] get wrapper from PyOP2

c3e84e1

spatial index for extruded meshes

5353b6c

initial support for extruded

78121f0

fall back on brute force for 1-dimension

ee9e28b

use fancy C++ kernel

3ef0e09

add libspatialindex dependency

3638bff

test case: point in cell query non-extruded

948ec73

test case: point in cell query extruded

7ac228e

test case: function spaces non-extruded

69d5b35

test case: function spaces extruded

185069b

[FFC] pass cdim to compile_element

4c44cc1

make cells affine for the function space test cases

d7c8d4e

n_cells -> n_cols

a423f75

refactor locate.cpp

94668f4

split spatialindex.pyx

c4f883d

spatialindex.pyx docstrings

b9ba392

refactor mesh.spatial_index

e3edd72

revise and test API

dc37d29

dorugeber reviewed Oct 16, 2015
View reviewed changes

miklos1 added 5 commits October 27, 2015 16:45

fix parallel evaluation bug

ecc520a

check input for parallel consistency

09ab80a

make parallel evaluation work (inefficiently)

1c0c35e

a parallel test

4804c3f

add section about point evaluation in the manual

88ab9f5

miklos1 force-pushed the point-evaluation-serial branch from 8f62f67 to 88ab9f5 Compare October 27, 2015 16:46

wence- reviewed Oct 28, 2015
View reviewed changes

miklos1 added 2 commits October 28, 2015 14:39

review fixes to code and documentation

e2d3d14

review fixes to test cases

6976c83

miklos1 force-pushed the point-evaluation-serial branch from 66b4fdd to 6976c83 Compare October 28, 2015 15:01

miklos1 merged commit 6976c83 into master Oct 28, 2015

miklos1 deleted the point-evaluation-serial branch October 28, 2015 17:04

This was referenced Oct 28, 2015

Point in cell queries #524

Closed

Evaluating a Function at a given point #155

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Point evaluations #612

Point evaluations #612

miklos1 commented Oct 16, 2015

dorugeber Oct 16, 2015

miklos1 commented Oct 20, 2015

dham commented Oct 20, 2015

miklos1 commented Oct 20, 2015

dham commented Oct 20, 2015

miklos1 commented Oct 21, 2015

miklos1 commented Oct 27, 2015

wence- Oct 28, 2015

wence- Oct 28, 2015

miklos1 Oct 28, 2015

wence- Oct 28, 2015

miklos1 Oct 28, 2015

miklos1 commented Oct 28, 2015

Point evaluations #612

Point evaluations #612

Conversation

miklos1 commented Oct 16, 2015

dorugeber Oct 16, 2015

Choose a reason for hiding this comment

miklos1 commented Oct 20, 2015

dham commented Oct 20, 2015

miklos1 commented Oct 20, 2015

dham commented Oct 20, 2015

miklos1 commented Oct 21, 2015

miklos1 commented Oct 27, 2015

wence- Oct 28, 2015

Choose a reason for hiding this comment

wence- Oct 28, 2015

Choose a reason for hiding this comment

miklos1 Oct 28, 2015

Choose a reason for hiding this comment

wence- Oct 28, 2015

Choose a reason for hiding this comment

miklos1 Oct 28, 2015

Choose a reason for hiding this comment

miklos1 commented Oct 28, 2015