Hello, Greetings!
I am getting error of Unable to allocate 75.2 GiB for an array with shape (100489, 100489) and data type float64. I am missing something fundamental in creating form with trial function dv. With this I am unable to create petsc.assemble_matrix of this bilinear form. Kindly help me to remove this high memory allocation error. Any help is highly appreciated!
The shell_curved_2elem.msh file used is:
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
6
1 -0.003241555532440543 0.002343569416552782 0.005201474763453007
2 -0.00235324539244175 0.003234538016840816 0.005203458480536938
3 -0.003245917148888111 0.00233752466738224 0.006653184536844492
4 -0.002358921570703387 0.00323040084913373 0.006654146127402782
5 -0.003808011068031192 0.001224357518367469 0.005199088249355555
6 -0.003809597576037049 0.001219412544742227 0.006647835951298475
$EndNodes
$Elements
2
1 3 2 2 2 1 2 4 3
2 3 2 2 2 3 6 5 1
$EndElements
The MWE code:
from dolfinx.io import gmshio
from mpi4py import MPI
import meshio
import dolfinx
from dolfinx.fem import VectorFunctionSpace, form, petsc, Function, FunctionSpace
from ufl import Jacobian, as_vector, sqrt, dot, cross, conditional, as_matrix
from ufl import lt,SpatialCoordinate, as_tensor, VectorElement, FiniteElement, MixedElement, Measure
from ufl import dx, dot, TrialFunction, TestFunction
import petsc4py.PETSc
mesh, subdomains, boundaries = gmshio.read_from_msh("shell_curved_2elem.msh", MPI.COMM_WORLD,0, gdim=3)
Info : Reading 'shell_curved_2elem.msh'...
Info : 6 nodes
Info : 2 elements
Info : Done reading 'shell_curved_2elem.msh'
#
# Creating local frame using linear shell tutorial
t = Jacobian(mesh)
t1 = as_vector([t[0, 0], t[1, 0], t[2, 0]])
t2 = as_vector([t[0, 1], t[1, 1], t[2, 1]])
e3 = cross(t1, t2)
e3 /= sqrt(dot(e3, e3))
ey = as_vector([0, 1, 0])
ez = as_vector([0, 0, 1])
e1 = cross(ey, e3)
norm_e1 = sqrt(dot(e1, e1))
e1 = conditional(lt(norm_e1, 0.5), ez, e1 / norm_e1)
e2 = cross(e3, e1)
norm_e2 = sqrt(dot(e2, e2))
e2 /= sqrt(dot(e2, e2))
# Creating ufl form
e=[e1,e2,e3]
def eps(v):
# taking derivative of v in local frame (local tangent gradient)
# dot(e1,as_vector([v[0].dx(i) for i in range(3)])) --> d(v[0])/d(e1)
ww=[dot(e[j],as_vector([v[0].dx(i) for i in range(3)])) for j in range(3)]
[w11,w12,w13]=ww
r21= dot(e2,as_vector(ww))
vv=e2[0]*r21
return vv
# creating form
V = VectorFunctionSpace(mesh, ("S", 1))
dv = TrialFunction(V)
v_ = TestFunction(V)
F=(eps(dv)*eps(v_))*dx
form(F)
---------------------------------------------------------------------------
MemoryError Traceback (most recent call last)
Cell In[5], line 39
36 v_ = TestFunction(V)
38 F=(eps(dv)*eps(v_))*dx
---> 39 form(F)
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:188, in form(form, dtype, form_compiler_options, jit_options)
185 return list(map(lambda sub_form: _create_form(sub_form), form))
186 return form
--> 188 return _create_form(form)
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:183, in form.<locals>._create_form(form)
180 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
181 return form argument"""
182 if isinstance(form, ufl.Form):
--> 183 return _form(form)
184 elif isinstance(form, collections.abc.Iterable):
185 return list(map(lambda sub_form: _create_form(sub_form), form))
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py:141, in form.<locals>._form(form)
139 if mesh is None:
140 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 141 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
142 form_compiler_options=form_compiler_options,
143 jit_options=jit_options)
145 # For each argument in form extract its function space
146 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
51 @functools.wraps(local_jit)
52 def mpi_jit(comm, *args, **kwargs):
53
54 # Just call JIT compiler when running in serial
55 if comm.size == 1:
---> 56 return local_jit(*args, **kwargs)
58 # Default status (0 == ok, 1 == fail)
59 status = 0
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
202 # Switch on type and compile, returning cffi object
203 if isinstance(ufl_object, ufl.Form):
--> 204 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.FiniteElementBase):
206 r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)
File /usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
197 except Exception:
198 pass
--> 199 raise e
201 obj, module = _load_objects(cache_dir, module_name, form_names)
202 return obj, module, (decl, impl)
File /usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
187 for name in form_names:
188 decl += form_template.format(name=name)
--> 190 impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
191 cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
192 except Exception as e:
193 try:
194 # remove c file so that it will not timeout next time
File /usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
256 import ffcx.compiler
258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
262 ffibuilder = cffi.FFI()
263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
264 extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)
File /usr/lib/python3/dist-packages/ffcx/compiler.py:102, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
100 # Stage 2: intermediate representation
101 cpu_time = time()
--> 102 ir = compute_ir(analysis, object_names, prefix, options, visualise)
103 _print_timing(2, time() - cpu_time)
105 # Stage 3: code generation
File /usr/lib/python3/dist-packages/ffcx/ir/representation.py:197, in compute_ir(analysis, object_names, prefix, options, visualise)
191 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
192 for e in analysis.unique_elements]
194 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
195 for e in analysis.unique_elements]
--> 197 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
198 options, visualise)
199 for (i, fd) in enumerate(analysis.form_data)]
200 ir_integrals = list(itertools.chain(*irs))
202 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
203 finite_element_names, dofmap_names, object_names)
204 for (i, fd) in enumerate(analysis.form_data)]
File /usr/lib/python3/dist-packages/ffcx/ir/representation.py:197, in <listcomp>(.0)
191 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
192 for e in analysis.unique_elements]
194 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
195 for e in analysis.unique_elements]
--> 197 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
198 options, visualise)
199 for (i, fd) in enumerate(analysis.form_data)]
200 ir_integrals = list(itertools.chain(*irs))
202 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
203 finite_element_names, dofmap_names, object_names)
204 for (i, fd) in enumerate(analysis.form_data)]
File /usr/lib/python3/dist-packages/ffcx/ir/representation.py:492, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise)
489 integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()}
491 # Build more specific intermediate representation
--> 492 integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
493 ir["entitytype"], integrands, ir["tensor_shape"],
494 options, visualise)
496 ir.update(integral_ir)
498 # Fetch name
File /usr/lib/python3/dist-packages/ffcx/ir/integral.py:84, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, p, visualise)
74 # Build terminal_data from V here before factorization. Then we
75 # can use it to derive table properties for all modified
76 # terminals, and then use that to rebuild the scalar graph more
77 # efficiently before argument factorization. We can build
78 # terminal_data again after factorization if that's necessary.
80 initial_terminals = {i: analyse_modified_terminal(v['expression'])
81 for i, v in S.nodes.items()
82 if is_modified_terminal(v['expression'])}
---> 84 mt_table_reference = build_optimized_tables(
85 quadrature_rule,
86 cell,
87 integral_type,
88 entitytype,
89 initial_terminals.values(),
90 ir["unique_tables"],
91 rtol=p["table_rtol"],
92 atol=p["table_atol"])
94 # Fetch unique tables for this quadrature rule
95 table_types = {v.name: v.ttype for v in mt_table_reference.values()}
File /usr/lib/python3/dist-packages/ffcx/ir/elementtables.py:366, in build_optimized_tables(quadrature_rule, cell, integral_type, entitytype, modified_terminals, existing_tables, rtol, atol)
364 # Clean up table
365 tbl = clamp_table_small_numbers(t['array'], rtol=rtol, atol=atol)
--> 366 tabletype = analyse_table_type(tbl)
368 if tabletype in piecewise_ttypes:
369 # Reduce table to dimension 1 along num_points axis in generated code
370 tbl = tbl[:, :, :1, :]
File /usr/lib/python3/dist-packages/ffcx/ir/elementtables.py:453, in analyse_table_type(table, rtol, atol)
450 elif is_ones_table(table, rtol=rtol, atol=atol):
451 # All values are 1.0
452 ttype = "ones"
--> 453 elif is_quadrature_table(table, rtol=rtol, atol=atol):
454 # Identity matrix mapping points to dofs (separately on each entity)
455 ttype = "quadrature"
456 else:
457 # Equal for all points on a given entity
File /usr/lib/python3/dist-packages/ffcx/ir/elementtables.py:420, in is_quadrature_table(table, rtol, atol)
418 def is_quadrature_table(table, rtol=default_rtol, atol=default_atol):
419 _, num_entities, num_points, num_dofs = table.shape
--> 420 Id = np.eye(num_points)
421 return (num_points == num_dofs and all(
422 np.allclose(table[0, i, :, :], Id, rtol=rtol, atol=atol) for i in range(num_entities)))
File /usr/lib/python3/dist-packages/numpy/lib/twodim_base.py:214, in eye(N, M, k, dtype, order, like)
212 if M is None:
213 M = N
--> 214 m = zeros((N, M), dtype=dtype, order=order)
215 if k >= M:
216 return m
MemoryError: Unable to allocate 75.2 GiB for an array with shape (100489, 100489) and data type float64