Hi dear Dooken,
Thank you for your help so far.
I just now want to study the problem in 1D and then visualize the solution using pyvista.
Following the discussions in https://fenicsproject.discourse.group/t/pyvista-with-mixed-elements-in-dolfinx/9082, I have done so:
import os
import matplotlib.pyplot as plt
import numpy as np
import ufl
from dolfinx.io import VTXWriter
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological, set_bc)
from dolfinx import log, plot
#from dolfinx.plot import vtk_mesh
from dolfinx.fem import Function, FunctionSpace
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, TrialFunction, derivative
from dolfinx import fem, io, mesh, plot
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from petsc4py import PETSc
try:
import pyvista as pv
import pyvistaqt as pvqt
have_pyvista = True
if pv.OFF_SCREEN:
pv.start_xvfb(wait=0.5)
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
# Save all logging to file
log.set_output_file("log.txt")
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
lx = 1.0 #length of x-axis
# parmeters
r_0 = 0.0
r_b = 1.0 # the coordinate of right boundary
msh = mesh.create_interval(MPI.COMM_WORLD, 10, points=(0, 1))
# Define ME function space
m1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
T1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
ME = FunctionSpace(msh, m1 * T1)
Q = FunctionSpace(msh, m1) #ME.sub(0)
V = FunctionSpace(msh, T1) #ME.sub(1)
# Defining our problem
am = fem.Constant(msh, ScalarType(3.e-06 ))
k = fem.Constant(msh, ScalarType(0.12))
c = fem.Constant(msh, ScalarType(1.)) #c = fem.Constant(msh, ScalarType(1284.))
rho_0 = fem.Constant(msh, ScalarType(1.)) #rho_0 = fem.Constant(msh, ScalarType(500.))
R = fem.Constant(msh, ScalarType(1.)) #R = fem.Constant(msh, ScalarType(2400.))
epsilon = fem.Constant(msh, ScalarType(0.1))
nu = fem.Constant(msh, ScalarType(1.e-02))
dt = 0.001 #5.0e-06 # 0.001 #0.02 #5.0e-06 # time step
#We create four python functions for determining the facets to apply boundary conditions to
def left(x):
return np.isclose(x[0], r_0) #left (x = 0)
def right(x):
return np.isclose(x[0], r_b) #right (x = 1)
boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_right)
bc_m_right = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_right, ME.sub(0))
boundary_facets_left = locate_entities_boundary(msh, msh.topology.dim - 1, left)
boundary_dofs_left = locate_dofs_topological(ME.sub(0), msh.topology.dim - 1, boundary_facets_left)
bc_m_left = dirichletbc(PETSc.ScalarType(0.12), boundary_dofs_left, ME.sub(0))
boundary_facets_right = locate_entities_boundary(msh, msh.topology.dim - 1, right)
boundary_dofs_right = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_right)
bc_T_right = dirichletbc(PETSc.ScalarType(30.), boundary_dofs_right, ME.sub(1))
boundary_facets_left = locate_entities_boundary(msh, msh.topology.dim - 1, left)
boundary_dofs_left = locate_dofs_topological(ME.sub(1), msh.topology.dim - 1, boundary_facets_left)
bc_T_left = dirichletbc(PETSc.ScalarType(30.), boundary_dofs_left, ME.sub(1))
#dirichlet boundary conditions
bcs = [bc_m_right, bc_m_left, bc_T_right, bc_T_left]
#Test Functions (q, v)
q, v = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u_n = Function(ME) # solution from previous converged step
#du = TrialFunction(ME)
# Split mixed functions
m, T = ufl.split(u)
m_n, T_n = ufl.split(u_n)
# Weak statement of the equations
m_n = Function(Q)
m_n.name = "m_n"
T_n = Function(V)
T_n.name = "T_n"
# The initial conditions are interpolated into a finite element space:
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
u.sub(0).interpolate(lambda x: np.full(x.shape[1], 0.25, dtype=np.float64)) #Initial condition for m
u.sub(1).interpolate(lambda x: np.full(x.shape[1], 10., dtype=np.float64)) #Initial condition for T
u.x.scatter_forward()
#Variational formulation of the system
F0 = inner(m, q) * dx - inner(m_n, q) * dx + dt*am * inner(grad(m), grad(q)) * dx + dt*am*nu * inner(grad(T), grad(q)) * dx
F1 = c*rho_0* inner(T, v) * dx - c*rho_0* inner(T_n, v) * dx + dt*k * inner(grad(T), grad(v)) * dx - R*epsilon * inner(m, v) * dx + R*epsilon * inner(m_n, v) * dx
F = F0 + F1
#Compute the jacobian
#Jac = derivative(F, u, du)
#Solve variational problem
# Create nonlinear problem and Newton solver
#problem = NonlinearProblem(F, u, bcs = bcs) #Replace
problem = NonlinearProblem(F, u, bcs = bcs)
#problem = NonlinearProblem(F, u, bcs = bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual" #residual #incremental
solver.rtol = 1e-12 #0.0 #1e-14 #0.0 #1e-12
#solver.max_it = 70
solver.report = True
#solver.max_it = 100 #To increase the number of iterations
solver.error_on_nonconvergence = False
solver.max_it = 1
info = solver.solve(u)
print(info)
with io.XDMFFile(MPI.COMM_WORLD, "problem_interval_mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
with dolfinx.io.VTXWriter(msh.comm, "moisture_solution1D.bp", [u.sub(0)], engine="BP4") as vtx:
vtx.write(0.0)
#with dolfinx.io.VTXWriter(msh.comm, "Temperature.bp", [u.sub(1).collapse()], engine="BP4") as vtx:
with dolfinx.io.VTXWriter(msh.comm, "Temperature_solution1D.bp", [u.sub(1)], engine="BP4") as vtx:
vtx.write(0.0)
import pyvista as pv
import matplotlib.pyplot as plt
pv.start_xvfb()
u0_h = u.sub(0).collapse() #m
u1_h = u.sub(1).collapse() #T
# Dummy function space for pyvista grid
W_dummy = FunctionSpace(msh, ("CG", 1))
#cells, types, x = plot.create_vtk_mesh(W_dummy) #dolfinx.plot.vtk_mesh
cells, types, x = dolfinx.plot.vtk_mesh(W_dummy) #dolfinx.plot.vtk_mesh
grid = pv.UnstructuredGrid(cells, types, x)
plotter = pv.Plotter()
plotter.open_gif("u0_htime.gif", fps=10)
grid.point_data["u0_h"] = u.sub(0).collapse().x.array.real
warped = grid.warp_by_scalar("u0_h", factor=1)
viridis = plt.colormaps.get_cmap("viridis").resampled(25)
sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
position_x=0.1, position_y=0.8, width=0.8, height=0.1)
renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
cmap=viridis, scalar_bar_args=sargs,
clim=[0, max(u.x.array)])
plotter.show_axes()
The code works but I don’t see any plot.
Please let me know what I am missing.
Thank you for your time and your patience.