Accessing entries of Global size matrix

General
1 
from fenics import *
import numpy as np

# Create mesh and function space
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'P', 1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
A = PETScMatrix()
assemble(a, tensor=A)

AM = np.array(A.array())
print(f"A[0,0] = {AM[0,0]}")

This is the basic idea that i am using to access few entries of the matrix A to update something for the next iteration in my code. But due to storing large size global matrix its giving me memory error. Can you provide me an alternative so that i don’t get this error and my code work for higher values of M like 256,384,512 etc.

2

Considering something along the lines of

from fenics import *
import numpy as np

# Create mesh and function space
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'P', 1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
A = PETScMatrix()
assemble(a, tensor=A)

def print_entries(entries):
    # DO NOT DO THIS IN AN ACTUAL CODE. You can use A.get instead
    AM = np.array(A.array())
    for entry in entries:
        print(f"Before: A[{entry[0], entry[1]}] = {AM[entry[0], entry[1]]}")

print("Before")
print_entries([(0, 0), (0, 1), (1, 0), (1, 1)])

A.set([[2.0, 3.0], [4.0, 5.0]], [0, 1], [0, 1])
A.apply("insert")

print("After")
print_entries([(0, 0), (0, 1), (1, 0), (1, 1)])

Furthermore read Read before posting: How do I get my question answered? (next time, please format the code properly) and The new DOLFINx solver is now recommended over DOLFIN

3

Thanks for your answer @francesco-ballarin. I want to update some values from this matrix to another matrix which is obtained by a loop over cells. The memory error is still ocurring due to use of A.array(). If i use A.get, then it is giving me the error TypeError: ‘method’ object is not subscriptable.

4

Slight modification of the above then

from fenics import *
import numpy as np

# Create mesh and function space
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'P', 1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
A = PETScMatrix()
assemble(a, tensor=A)

def print_entries_inefficient(entries):
    # DO NOT DO THIS IN AN ACTUAL CODE. You can use A.get instead
    AM = np.array(A.array())
    for entry in entries:
        print(f"Before: A[{entry[0], entry[1]}] = {AM[entry[0], entry[1]]}")

def print_entries(entries0, entries1):
    values = np.zeros((len(entries0), len(entries1)))
    A.get(values, entries0, entries1)
    for idx0, entry0 in enumerate(entries0):
        for idx1, entry1 in enumerate(entries1):
            print(f"Before: A[{entry0, entry1}] = {values[idx0, idx1]}")

print("Before (inefficient)")
print_entries_inefficient([(0, 0), (0, 1), (1, 0), (1, 1)])
print("Before")
print_entries([0, 1], [0, 1])

A.set([[2.0, 3.0], [4.0, 5.0]], [0, 1], [0, 1])
A.apply("insert")

print("After (inefficient)")
print_entries_inefficient([(0, 0), (0, 1), (1, 0), (1, 1)])
print("After")
print_entries([0, 1], [0, 1])
1 Like