Git Product home page Git Product logo

torch-fenics's Introduction

Torch-FEniCS

The torch-fenics package enables models defined in FEniCS to be used as modules in PyTorch.

Install

Install FEniCS and run

pip install git+https://github.com/barkm/torch-fenics.git@master

A clean install of the package and its dependencies can for example be done with Conda

conda create --name torch-fenics
conda activate torch-fenics
conda install -c conda-forge fenics
pip install git+https://github.com/barkm/torch-fenics.git@master

Details

FEniCS objects are represented in PyTorch using their corresponding vector representation. For finite element functions this corresponds to their coefficient representation.

The package relies on dolfin-adjoint in order for the FEniCS module to be compatible with the automatic differentiation framework in PyTorch

Example

The torch-fenics package can for example be used to define a PyTorch module which solves the Poisson equation using FEniCS.

The process of solving the Poisson equation in FEniCS can be specified as a PyTorch module by deriving the torch_fenics.FEniCSModule class

# Import fenics and override necessary data structures with fenics_adjoint
from fenics import *
from fenics_adjoint import *

import torch_fenics

# Declare the FEniCS model corresponding to solving the Poisson equation
# with variable source term and boundary value
class Poisson(torch_fenics.FEniCSModule):
    # Construct variables which can be in the constructor
    def __init__(self):
        # Call super constructor
        super().__init__()

        # Create function space
        mesh = UnitIntervalMesh(20)
        self.V = FunctionSpace(mesh, 'P', 1)

        # Create trial and test functions
        u = TrialFunction(self.V)
        self.v = TestFunction(self.V)

        # Construct bilinear form
        self.a = inner(grad(u), grad(self.v)) * dx

    def solve(self, f, g):
        # Construct linear form
        L = f * self.v * dx

        # Construct boundary condition
        bc = DirichletBC(self.V, g, 'on_boundary')

        # Solve the Poisson equation
        u = Function(self.V)
        solve(self.a == L, u, bc)

        # Return the solution
        return u

    def input_templates(self):
        # Declare templates for the inputs to Poisson.solve
        return Constant(0), Constant(0)

The Poisson.solve function can now be executed by giving the module the appropriate vector input corresponding to the input templates declared in Poisson.input_templates. In this case the vector representation of the template Constant(0) is simply a scalar.

# Construct the FEniCS model
poisson = Poisson()

# Create N sets of input
N = 10
f = torch.rand(N, 1, requires_grad=True, dtype=torch.float64)
g = torch.rand(N, 1, requires_grad=True, dtype=torch.float64)

# Solve the Poisson equation N times
u = poisson(f, g)

The output of the can now be used to construct some functional. Consider summing up the coefficients of the solutions to the Poisson equation

# Construct functional 
J = u.sum()

The derivative of this functional with respect to f and g can now be computed using the torch.autograd framework.

# Execute backward pass
J.backward() 

# Extract gradients
dJdf = f.grad
dJdg = g.grad

Developing

Install dependencies

conda env create -n torch-fenics -f environment.yml
conda activate torch-fenics

Install package in editable mode

pip install -e .[test]

The unit-tests can then be run as follows

python -m pytest tests

torch-fenics's People

Contributors

mossjacob avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

torch-fenics's Issues

Issue with the poisson example

Hi !
I tried to install your package and to test your poisson example, but an issue occured.
This is how I installed everything:

 conda create --name fenics_torch
 source activate fenics_torch
 conda install -c conda-forge fenics=2018.1.0
 pip install --process-dependency-links git+https://github.com/pbarkm/torch-fenics.git@master

And this is the output of python poisson.py:

/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/ffc/uflacs/analysis/dependencies.py:61: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  active[targets] = 1
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/ffc/uflacs/analysis/dependencies.py:61: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  active[targets] = 1
Traceback (most recent call last):
  File "poisson.py", line 69, in <module>
    J.backward()
  File "/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/torch/tensor.py", line 93, in backward
    torch.autograd.backward(self, gradient, retain_graph, create_graph)
  File "/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/torch/autograd/__init__.py", line 90, in backward
    allow_unreachable=True)  # allow_unreachable flag
  File "/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/torch/autograd/function.py", line 76, in apply
    return self._forward_cls.backward(self, *args)
  File "/Users/user/anaconda2/envs/fenics_torch/lib/python3.6/site-packages/torch_fenics/torch_fenics.py", line 127, in backward
    adj_value=adj_value)
TypeError: compute_gradient() got an unexpected keyword argument 'adj_value'

AttributeError: 'Function' object has no attribute 'block_variable'

I have written a subclass of torch_fenics. In this, the input is a vector from DG space. I use this input in the weak formulation and then calculate the solution. Further, I need the gradient of the solution with respect to the given input. The code to reproduce is given below followed by the error log. Please help me in fixing this issue.

import torch
from fenics import *
import torch_fenics
from fenics_adjoint import *
from dolfin import *

class Test(torch_fenics.FEniCSModule):
   
    def __init__(self,nx=5,ny=5): 
        super(Test,self).__init__()
        mesh = UnitSquareMesh(nx,ny)
        self.V = FunctionSpace(mesh, 'P', 1)
        self.W = FunctionSpace(mesh, 'DG', 0)
        self.u = TrialFunction(self.V)
        self.v = TestFunction(self.V)
        self.b = Constant((2.0, 3.0))
        
    def solve(self,inp):
        
        def E0_boundary(x, on_boundary):
            return on_boundary

        bc = DirichletBC(self.V, Constant(0.0),E0_boundary)
        u_sol = Function(self.V)
        a = (inp*(dot(self.b,grad(self.u)))*dot(self.b,grad(self.v)))*dx()  
        L = (inp*dot(self.b,grad(self.v)))*dx()
        solve(a == L, u_sol, bcs = bc)
        return u_sol
        
    def input_templates(self):
        return Function(self.W)
        
nx,ny = 5,5
n_cells = nx*ny*2
inp = torch.rand(1,n_cells,requires_grad = True,dtype = torch.float64)
test = Test(nx,ny)
usol = test(inp)
loss = torch.norm(usol)
loss.backward()

I get the following error log on running the same

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-1e15af964879> in <module>
     37 usol = test(tau)
     38 loss = torch.norm(usol)
---> 39 loss.backward()

~/miniconda3/envs/py37/lib/python3.7/site-packages/torch/tensor.py in backward(self, gradient, retain_graph, create_graph)
    193                 products. Defaults to ``False``.
    194         """
--> 195         torch.autograd.backward(self, gradient, retain_graph, create_graph)
    196 
    197     def register_hook(self, hook):

~/miniconda3/envs/py37/lib/python3.7/site-packages/torch/autograd/__init__.py in backward(tensors, grad_tensors, retain_graph, create_graph, grad_variables)
     97     Variable._execution_engine.run_backward(
     98         tensors, grad_tensors, retain_graph, create_graph,
---> 99         allow_unreachable=True)  # allow_unreachable flag
    100 
    101 

~/miniconda3/envs/py37/lib/python3.7/site-packages/torch/autograd/function.py in apply(self, *args)
     75 
     76     def apply(self, *args):
---> 77         return self._forward_cls.backward(self, *args)
     78 
     79 

~/miniconda3/envs/py37/lib/python3.7/site-packages/torch_fenics/torch_fenics.py in backward(ctx, *grad_outputs)
     88         # Check which gradients need to be computed
     89         controls = list(map(fenics_adjoint.Control,
---> 90                             (c for g, c in zip(ctx.needs_input_grad[1:], ctx.fenics_inputs) if g)))
     91 
     92         # Compute and accumulate gradient for each output with respect to each input

~/miniconda3/envs/py37/lib/python3.7/site-packages/pyadjoint/control.py in __init__(self, control)
     38     def __init__(self, control):
     39         self.control = control
---> 40         self.block_variable = control.block_variable
     41 
     42     def data(self):

AttributeError: 'Function' object has no attribute 'block_variable'

Optimizing boundary condition for given 2D solution

Hello,
I am trying to use torch-fenics to compute gradients with respect to input array which is imposed as a Dirichlet boundary condition on the top surface of the domain. I have written a simple solver in FEniCS which extrudes the given 1 dimensional input to a 2D domain by solving du/dy = 0, subject to boundary condition u(x) = some input array on the top boundary. To achieve this, I create a 1D function space using SubMesh on top boundary of the 2D mesh and I use this 1D function space to represent the input function in input_templates() for torch_fenics module. Then I pass a 1D torch tensor as input to the solver. Here is a minimal example of my implementation:

import numpy as np
import dolfin as dl
import torch
from fenics import *
from fenics_adjoint import *
import torch_fenics

class Top(SubDomain):
	def inside(self, x, on_boundary):
		return near(x[1], 1.0, DOLFIN_EPS)

class Extrude(torch_fenics.FEniCSModule):
    def __init__(self):
        super().__init__()
		# define the mesh and mesh function
        self.mesh = UnitSquareMesh(nx=10, ny=10)
        self.mesh_function = MeshFunction("size_t", self.mesh, 1)
        self.mesh_function.set_all(0)
        
		# mark top boundary
        self.top = Top()
        self.top.mark(self.mesh_function, 1)
        
        self.scalar_function_space = FunctionSpace(self.mesh, "Lagrange", 1)

		# create 1D function space corresponding to top boundary top used for input_templates
        bm = BoundaryMesh(self.mesh, "exterior")
        self.mesh_1d = SubMesh(bm, self.top)
        self.function_space_1d = FunctionSpace(self.mesh_1d, "Lagrange", 1)
        self.u = TrialFunction(self.scalar_function_space)
        self.v = TestFunction(self.scalar_function_space)
        self.u_sol = Function(self.scalar_function_space)
    
    def solve(self, T_in):
        T_in.set_allow_extrapolation(True)
        self.T_in = T_in
        self.bc1 = DirichletBC(self.scalar_function_space, self.T_in, self.top, "pointwise")

        # define the weak for du/dy = 0
        weak_form = grad(self.u)[1] * self.v * dx
        weak_lhs = lhs(weak_form)
        weak_rhs = rhs(weak_form)

        solve(weak_lhs == weak_rhs, self.u_sol, [self.bc1])
        
        return self.u_sol

    def input_templates(self):
        return Function(self.function_space_1d)

extrude = Extrude()

# set the input array (top boundary values)
u_1d = np.linspace(1.0, 0.0, 11, dtype=np.float64) + 1.0
u_1d = torch.tensor(u_1d, requires_grad=True, dtype=torch.float64).reshape(1,-1)

u_2d = extrude(u_1d)
J = u_2d.sum()
print(f"J = {J}")
J.backward()
print(u_1d.grad)

The solver seems to work fine and it extrudes the input array to the 2D domain as expected, but when calling J.backward(), I get the following error saying FunctionSpaces are expected to be on same mesh:

J = 181.5
Traceback (most recent call last):
  File "test.py", line 60, in <module>
    J.backward()
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/torch/_tensor.py", line 396, in backward
    torch.autograd.backward(self, gradient, retain_graph, create_graph, inputs=inputs)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/torch/autograd/__init__.py", line 175, in backward
    allow_unreachable=True, accumulate_grad=True)  # Calls into the C++ engine to run the backward pass
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/torch/autograd/function.py", line 253, in apply
    return user_fn(self, *args)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/torch_fenics/torch_fenics.py", line 96, in backward
    tape=ctx.tape, adj_value=adj_value)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/pyadjoint/drivers.py", line 28, in compute_gradient
    tape.evaluate_adj(markings=True)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/pyadjoint/tape.py", line 140, in evaluate_adj
    self._blocks[i].evaluate_adj(markings=markings)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/pyadjoint/tape.py", line 47, in wrapper
    return function(*args, **kwargs)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/pyadjoint/block.py", line 134, in evaluate_adj
    prepared)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/fenics_adjoint/types/dirichletbc.py", line 143, in evaluate_adj_component
    r = compat.extract_bc_subvector(adj_value, c.function_space(), bc)
  File "~/anaconda3/envs/fenics/lib/python3.7/site-packages/fenics_adjoint/types/compat.py", line 218, in extract_bc_subvector
    assigner = backend.FunctionAssigner(Vtarget, backend.FunctionSpace(bc.function_space()))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [email protected]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create function assigner.
*** Reason:  Expected all FunctionSpaces to be defined over the same Mesh.
*** Where:   This error was encountered inside FunctionAssigner.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  c5b9b269f4a6455a739109e3a66e036b5b8412f5
*** -------------------------------------------------------------------------

I am quite new to both FEniCS and Torch-FEniCS, and I would greatly appreciate any help in solving this issue. Also, please let me know if there is a better way to achieve this kind of optimization with respect to the boundary condition than creating 1D SubMesh and FunctionSpace as input.

Thanks in advance,
Shashwat

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.