February 2022 - firedrake - imperial.ac.uk

Error with Point Evaluation and VertexOnlyMesh
by Venketeswaran, Abhishek (FN) 01 Mar '22

01 Mar '22

Dear Koki and Lawrence, I am having trouble extracting my results at certain points on my mesh. I understand that there are 2 methods to do this (i) at() function and (ii) VertexOnlyMesh. Both dont work for me. While the former raises Point Not in Domain Error, the latter raises TypeError: an integer is required. Here is a minimal example of trying both methods out: from firedrake import * from firedrake.petsc import PETSc from firedrake.mesh import VertexOnlyMesh from collections import OrderedDict import numpy as np mesh = UnitSquareMesh(10, 10) V = FunctionSpace(mesh, "CG", 1) v = TestFunction(V) u = TrialFunction(V) a = inner(grad(u), grad(v)) * dx L = Constant(0) * v * dx bc1 = DirichletBC(V, 0.0, 1) bc2 = DirichletBC(V, 1.0, 2) sol = Function(V) problem1 = LinearVariationalProblem(a, L, sol, bcs=[bc1, bc2], constant_jacobian=True) solver1 = LinearVariationalSolver(problem1, solver_parameters={"ksp_type": "cg", "pc_type": "jacobi"}) solver1.solve() print(sol.dat.data_ro) crd=[0.5256,0.5051] num_points = 25 δs = np.linspace(0.0, 1.0, num_points + 1) X, Y = np.meshgrid(δs, δs) xs = np.vstack((X.flatten(), Y.flatten())).T vm=1 pe=1 #Test VertexOnlyMesh Generation if vm: point_cloud = VertexOnlyMesh(mesh, xs) #Test Point Evaluation of Function if pe: print(sol.at(crd)) Kindly let me know where I should start looking to fix this problem. -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran(a)netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:b3222045-e3c2-4f40-9a6f-b6c662401b72]

1 1

Obtaining the solution for a subset of nodes in Parallel
by Venketeswaran, Abhishek (FN) 12 Feb '22

12 Feb '22

Hi, I am running a simulation on a mesh with large number of nodes. I need to capture the results at a subset of the nodes. I am running into trouble when I execute the job in parallel. Kindly take a look at the lines in red to understand the problem. While I am not getting an error, I am also not able to get the correct results here. I would appreciate your assistance in this regard. Thanks Abhishek ########################################### # Obtaining the solution for a subset of mesh nodes ########################################## x,y,z = SpatialCoordinate(W.mesh()) phi=ufl.mathfunctions.Atan2(y,x) Nsteps=int(T/deltat) time=np.linspace(0,T,Nsteps+1) resnodesx,rescoords=axialslice(mesh,Ro,np.pi/2) # axialslice returns a set of nodes and their coordinates corresponding to a particular criterion resnodes=resnodesx.astype(np.int32) #PETSc gives an error if the indices are not int32 # Trying to save the solution components at all time instants for a particular subset uradial=np.zeros((len(resnodes),len(time)+3)) utang=np.zeros((len(resnodes),len(time)+3)) uaxial=np.zeros((len(resnodes),len(time)+3)) uradial[:,:3]=rescoords utang[:,:3]=rescoords uaxial[:,:3]=rescoords for i in range(len(time)): t = time[i] collar_disp.assign([0.0,0.0,collar_load(t,T_pulse,f_c,amp=dispamp)]) solver1.solve() # Update old fields with new quantities u_nm1.assign(u_n) u_n.assign(u_np1) u_r.assign(interpolate(u_n[0]*cos(phi)+u_n[1]*sin(phi),W)) u_t.assign(interpolate(-u_n[0]*sin(phi)+u_n[1]*cos(phi),W)) u_z.assign(interpolate(u_n[2],W)) with u_r.dat.vec_ro as vur: #Obtained this idea from Printing in Parallel page of Firedrake Docs urad=vur.getValues(resnodes) print(len(urad)) # length of urad seems to differ for each process uradial[:,i+3]=urad with u_t.dat.vec_ro as vut: utan=vut.getValues(resnodes) utang[:,i+3]=utan with u_z.dat.vec_ro as vuz: uax=vuz.getValues(resnodes) uaxial[:,i+3]=uax if i%20 ==0: outfile.write(u_n,u_r,u_t,u_z, time=t) gc.collect() #Consequently the CSV's contain only a subset of the nodal results np.savetxt(csvradfilename,uradial,fmt='%1.6e',delimiter=',') np.savetxt(csvtangfilename,utang,fmt='%1.6e',delimiter=',') np.savetxt(csvaxialfilename,uaxial,fmt='%1.6e',delimiter=',') -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran(a)netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:fc8c5c43-be1a-42c1-880c-cb38951ae7af]

1 0

Firedrake support on GitHub Discussions
by Ham, David A 02 Feb '22

02 Feb '22

Dear Firedrakers, I’m writing to update you on a change to the preferred support routes for Firedrake. The preferred way to ask for Firedrake help is now on GitHub Discussions. This will enable the Firedrake community to build up a persistent and searchable collection of answers to Firedrake questions, and means that users only need a GitHub account to post questions. Please post your requests here: https://github.com/firedrakeproject/firedrake/discussions/categories/firedr… Please also always feel welcome to discuss your issues with us on the Firedrake Slack channel. We always enjoy interacting with and helping our users. Yours, David -- Reader in Computational Mathematics Department of Mathematics Imperial College London [Graphical user interface, text Description automatically generated]<https://www.amazon.co.uk/dp/B09NR9R9P3>

1 0

Help with implementing steps in simulation
by Venketeswaran, Abhishek (FN) 01 Feb '22

01 Feb '22

Hi, I would like to prescribe a boundary condition for an initial time duration and then remove it for the remaining time duration. Commercial finite element solvers allow you to do this by specifying steps. I tried to implement it using the code below. This is giving nonsensical results at and after the time at which the BCs change. Can you please assist me in how I may implement this properly. I can share the full code if needed. m = rho*inner((u - 2.0 * u_n + u_nm1),v) / Constant(dt * dt) * dxlump k= inner(sigma(u_n),epsilon(v))*dx def collar_load(t,T,freq,amp=1.0): return amp*hanning(t,T)*stepfn(t,T)*np.sin(freq*t) # Define the loading as an expression depending on t t=0 collar_disp= Constant([0.0, 0.0, 0.0]) F = m+k a, r = lhs(F), rhs(F) zero = Constant((0.0, 0.0, 0.0)) # Dirichlet BC that is valid for the entire time duration bc1 = DirichletBC(V, zero, 11) # Collar BC that remains on till T_pulse bc2= DirichletBC(V, collar_disp, 21) # My attempt of defining two problems to mimic two simulation steps problem1 = LinearVariationalProblem(a, r, u_np1, bcs=[bc1,bc2], constant_jacobian=True) solver1 = LinearVariationalSolver(problem1, solver_parameters=params) problem2 = LinearVariationalProblem(a, r, u_np1, bcs=[bc1], constant_jacobian=True) solver2 = LinearVariationalSolver(problem2, solver_parameters=params) for i in range(len(time)): t = time[i] if (t<=T_pulse): solver1.solve() else: solver2.solve() u_nm1.assign(u_n) u_n.assign(u_np1) -- Abhishek Venketeswaran NETL Research Associate - ORISE National Energy Technology Laboratory Department of Energy abhishek.venketeswaran(a)netl.doe.gov Work: 412-386-4833 Mobile: 716-507-7890 www.netl.doe.gov<http://www.netl.doe.gov> [cid:97ebd409-8c9c-4361-9753-e4c626433151]

2 1