From 23431e78e45429a9db1973fee3c13586e13d4c49 Mon Sep 17 00:00:00 2001 From: Reuben Hill Date: Fri, 14 Feb 2020 16:38:23 +0000 Subject: [PATCH 1/3] first assembly demo commit --- assembly_demo.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 assembly_demo.py diff --git a/assembly_demo.py b/assembly_demo.py new file mode 100644 index 0000000..68b186e --- /dev/null +++ b/assembly_demo.py @@ -0,0 +1,38 @@ +import numpy as np +from firedrake import * + +# set up a really dull mesh and function space +m = UnitSquareMesh(5, 5) +V = FunctionSpace(m, 'CG', 1) +# Define a Function in the FunctionSpace - this is a UFL `Coefficient` +f = Function(V) +# Get symbolic representation of coordinates in mesh +x, y = SpatialCoordinate(m) +# Set f to be some really boring function (here a plane) +f.interpolate(2.*x) +# Define an integral of the boring function over the interior of the +# mesh as a UFL form. This is a 0-form since there are no unknown +# functions (UFL `Argument`s) only the known function f (a UFL +# `Coefficient`) +zeroform = f*dx +# Call assemble to evaluate this 0-form as a rank zero tensor (a real +# number). +# This is done in two steps: +# - 1. assemble._assemble is called to yield callable functions +# (typically PYOP2 loops) +# - a. tsfc_interface.compile_form is used to generate kernels to +# execute on the mesh with PYOP2 +# - b. PYOP2 is configured to execute these loops, with the +# functions to execute to make this happen being yielded by +# assemble._assemble +# - 2. the loops are called to yield the final result +assembled = assemble(zeroform) +assert np.isclose(assembled, 1.) + +# Now let's add an unknown TestFunction (a special case of a UFL +# `Argument`) to our form to create a 1-form +v = TestFunction(V) +oneform = v*dx +# Call assemble to evaluate this 1-form as a rank one tensor (a vector) +assembled = assemble(oneform) +print(assembled) \ No newline at end of file From 7b885b3a9c159939a9f26244b4195af57b5892a6 Mon Sep 17 00:00:00 2001 From: Reuben Hill Date: Fri, 20 Mar 2020 15:30:41 +0000 Subject: [PATCH 2/3] update with CJC input --- assembly_demo.py | 58 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 49 insertions(+), 9 deletions(-) diff --git a/assembly_demo.py b/assembly_demo.py index 68b186e..40165d7 100644 --- a/assembly_demo.py +++ b/assembly_demo.py @@ -1,6 +1,10 @@ import numpy as np +from matplotlib import pyplot as plt from firedrake import * +# In the simplest case you can use assemble to evaluate integrals with no +# unknowns (0-forms): + # set up a really dull mesh and function space m = UnitSquareMesh(5, 5) V = FunctionSpace(m, 'CG', 1) @@ -13,26 +17,62 @@ # Define an integral of the boring function over the interior of the # mesh as a UFL form. This is a 0-form since there are no unknown # functions (UFL `Argument`s) only the known function f (a UFL -# `Coefficient`) +# `Coefficient`). zeroform = f*dx # Call assemble to evaluate this 0-form as a rank zero tensor (a real # number). # This is done in two steps: # - 1. assemble._assemble is called to yield callable functions -# (typically PYOP2 loops) +# (typically PyOP2 loops) # - a. tsfc_interface.compile_form is used to generate kernels to -# execute on the mesh with PYOP2 -# - b. PYOP2 is configured to execute these loops, with the +# execute on the mesh with PyOP2 +# - b. PyOP2 is configured to execute these loops, with the # functions to execute to make this happen being yielded by # assemble._assemble -# - 2. the loops are called to yield the final result +# - 2. the loops are called to yield the final result - here a real number. assembled = assemble(zeroform) assert np.isclose(assembled, 1.) -# Now let's add an unknown TestFunction (a special case of a UFL -# `Argument`) to our form to create a 1-form +# Now consider a linear variational problem within a finite element +# function space V: find u in V such that a(u,v) = F(v) where v is in V. +# Given that `u = \sum_i u_i \phi_i` and `v = sum_i v_i\phi_i` where phi_i +# are the n basis functions in V, we can turn this into an equivalent matrix- +# vector equation +# +# `\sum_j a(\phi_i, \phi_j)u_j = F(\phi_i) i = 1, ..., n` +# +# where \phi_i are known basis functions and u_j are the unknown +# coefficients of u. + +# Lets look at the linear form F(v) first. +# Note that `TestFunction` constructs a special case of a UFL `Argument` +# defined on our function space. v = TestFunction(V) +# Since v is unknown, F(v) is a 1-form: a linear map from the test +# function v to a set of n numbers `F(\phi_i)` where i = 1, ..., n +# (a rank 1 tensor). +# Below is the simplest possible example: oneform = v*dx -# Call assemble to evaluate this 1-form as a rank one tensor (a vector) +# Assemble works as for the 0-form but with the final result being the +# rank 1 tensor `F(\phi_i)` where i = 1, ..., n. assembled = assemble(oneform) -print(assembled) \ No newline at end of file +# Note that in Firedrake, we reuse the `Function` type to store each +# `F(\phi_i)` where we usually store the basis coefficient u_i +# corresponding to basis function `\phi_i`. +print(type(assembled)) + +# Next lets look at the bilinear form a(u,v). +# `TrialFunction` also constructs a special case of a UFL `Argument` +# defined on our function space. +u = TrialFunction(V) +# Since u and v are unknown a(u,v) is a 2-form: a bilinear map from the +# trial function u and test function v to a set of numbers `a(\phi_i, \phi_j)` +# where i = 1, ..., n and j = 1, ..., n (a rank 2 tensor). +# Below is the simplest possible example: +twoform = v*u*dx +# Assemble again works as for the 0-form but with the final result being +# the rank 2 tensor `a(\phi_i, \phi_j)` where i = 1, ..., n and j = 1, ..., n +assembled = assemble(twoform) +# Assembled matrices are saved in their own `Matrix` type. +print(type(assembled)) + From 0a9fb1168f15a50c99ae5ad1094f6a189c2b66d0 Mon Sep 17 00:00:00 2001 From: Reuben Hill Date: Fri, 20 Mar 2020 15:51:08 +0000 Subject: [PATCH 3/3] add solve --- assembly_demo.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/assembly_demo.py b/assembly_demo.py index 40165d7..7316c7f 100644 --- a/assembly_demo.py +++ b/assembly_demo.py @@ -30,8 +30,8 @@ # functions to execute to make this happen being yielded by # assemble._assemble # - 2. the loops are called to yield the final result - here a real number. -assembled = assemble(zeroform) -assert np.isclose(assembled, 1.) +number = assemble(zeroform) +assert np.isclose(number, 1.) # Now consider a linear variational problem within a finite element # function space V: find u in V such that a(u,v) = F(v) where v is in V. @@ -55,11 +55,11 @@ oneform = v*dx # Assemble works as for the 0-form but with the final result being the # rank 1 tensor `F(\phi_i)` where i = 1, ..., n. -assembled = assemble(oneform) +assembled_rhs_vector = assemble(oneform) # Note that in Firedrake, we reuse the `Function` type to store each # `F(\phi_i)` where we usually store the basis coefficient u_i # corresponding to basis function `\phi_i`. -print(type(assembled)) +print(type(assembled_rhs_vector)) # Next lets look at the bilinear form a(u,v). # `TrialFunction` also constructs a special case of a UFL `Argument` @@ -72,7 +72,19 @@ twoform = v*u*dx # Assemble again works as for the 0-form but with the final result being # the rank 2 tensor `a(\phi_i, \phi_j)` where i = 1, ..., n and j = 1, ..., n -assembled = assemble(twoform) +assembled_lhs_matrix = assemble(twoform) # Assembled matrices are saved in their own `Matrix` type. -print(type(assembled)) +print(type(assembled_lhs_matrix)) +# let's solve the linear system, placing our solution into a newly made +# `Function` x +A = assembled_lhs_matrix +x = Function(V) +b = assembled_rhs_vector +solver = LinearSolver(A) +solver.solve(x, b) +u = x +# Note that the linear variation problem we solved has solution u = 1. +# Let's see how good our solution is: +plot(u) +plt.show() \ No newline at end of file