Hello,
I’m trying to simulate the melting of a glacier.
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
# Define the parameters of the problem
Lx = 1000.0 # Length of the domain in x direction (m)
Ly = 1000.0 # Length of the domain in y direction (m)
Lz = 500.0 # Length of the domain in z direction (m)
nx = 100 # Number of elements in x direction
ny = 100 # Number of elements in y direction
nz = 50 # Number of elements in z direction
alpha = 1.0e-6 # Thermal diffusivity (m^2/s)
L = 3.35e5 # Latent heat of fusion (J/kg)
rho = 917.0 # Density of ice (kg/m^3)
Cp = 2060.0 # Specific heat capacity of ice (J/kg-K)
Tm = 273.15 # Melting point of ice (K)
T0 = 253.15 # Initial temperature (K)
t_max = 5.0 * 365.0 * 24.0 * 3600.0 # Maximum simulation time (s)
dt = 10000.0 # Time step (s)
# Define the mesh
mesh = BoxMesh(Point(0, 0, 0), Point(Lx, Ly, Lz), nx, ny, nz)
# Define the function spaces
P1 = FiniteElement('P', tetrahedron, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
# Define the test functions
v, q = TestFunctions(V)
# Define the functions for the temperature and phase fields
u = Function(V)
u0 = Function(V)
# Split the mixed functions
T, phi = split(u)
T0, phi0 = split(u0)
# Define the boundary conditions
def boundary(x, on_boundary):
return on_boundary
bc_T = DirichletBC(V.sub(0), T0, boundary) #1 https://fenicsproject.org/qa/11979/significance-of-collapsing/
bc_phi = DirichletBC(V.sub(1), phi0, boundary)
bc = [bc_T, bc_phi]
# Define the initial conditions
T.interpolate(Constant(T0))
phi.interpolate(Constant(0.0))
phi.vector()[T.vector() > Tm] = 1.0
# Define the weak forms for the temperature and phase fields
a_T = (1/dt)*inner(T - T0, v)*dx + alpha*inner(grad(T), grad(v))*dx
L_T = (1/dt)*inner(T0, v)*dx
a_phi = (1/dt)*inner(phi - phi0, q)*dx + inner(grad(phi), grad(q))*dx + M(phi)/epsilon*inner(grad(T), grad(q))*dx
L_phi = (1/dt)*inner(phi0, q)*dx
# Define the time-stepping loop
t = 0.0
while t < t_max:
# Solve the coupled system of equations
solve(a_T == L_T, T, bc_T)
solve(a_phi == L_phi, phi, bc_phi)
# Update the temperature and phase fields
T0.assign(T)
phi0.assign(phi)
# Update the time
t += dt
# Plot the results
if int(t/dt) % 100 == 0:
plt.figure()
plt.subplot
same code but with line numbering:
But I came across this error message while executing the code.
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
Traceback (most recent call last):
File "/home/edward/Documents/Projects/python code with finite element method for Melting of a glacier example/main.py", line 43, in <module>
bc_T = DirichletBC(V.sub(0), T0, boundary) #1 https://fenicsproject.org/qa/11979/significance-of-collapsing/
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/dirichletbc.py", line 76, in __init__
expr = project(args[1], args[0]) # FIXME: This should really be interpolaton (project is expensive)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 137, in project
function = Function(V)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py", line 220, in __init__
self._cpp_object = cpp.function.Function(V._cpp_object)
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
***
*** fenics-support@googlegroups.com
***
*** 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.
*** Reason: Cannot be created from subspace. Consider collapsing the function space.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
I have little to no experience with Fenics and would appreciate any help.