Source term as the curl of another field? Second order transient term?

[yurigba]

Hi,

I am trying to write something similar to a leap-frog scheme in FiPy. I am trying to understand two things:

1. I want to implement the curl of another vector field as a source term. How this should be implemented in FiPy?

The values for the curl comes from the cell vertexes. The curl should apply to a rank-1 tensor for now, and the first problem is that in 2D problems the curl needs to be adapted properly such that if the field is orthogonal to the 2D plane the result should be on the 2D plane, and if the field is in the 2D plane, the curl should be orthogonal to the plane.

If there is not a way to do this using a term module, is there a way to directly input the values from the other CellVariable?

2. Is there anything like a second order time derivative transient term? If there is not, how is the "right way" of creating an update formula that is like the second order time derivative?

The values from the old step should come from cell center but also from adjacent cells. This will be needed for a higher-order scheme that I need to implement. If the needed modules are not implemented, I plan on making my own transient terms and curl terms (located on face-values), and if it goes well with FiPy I would like to contribute with them.

Thanks in advance!

[wd15]

1. I want to implement the curl of another vector field as a source term. How this should be implemented in FiPy?

I'm not sure why you'd want this as a source term. Isn't the variable being operated on also the solution variable?

If there is not a way to do this using a term module, is there a way to directly input the values from the other CellVariable?

I think that having this has a term makes the most sense.

2. Is there anything like a second order time derivative transient term? If there is not, how is the "right way" of creating an update formula that is like the second order time derivative?

You can split the second order time derivative into two equations like for the wave equation. This makes the most sense I think.

Are you interested in purely hyperbolic systems?

[yurigba]

Are you interested in purely hyperbolic systems?

Yes. Specifically I want to do a finite-volume formulation of Maxwell's equations, and FiPy has the tools to make this less cumbersome, if I learn how to use it properly. So there are two coupled equations, Faraday's Law and Ampère's Law. So the way I understand that would work in FiPy is to work the time derivative of E as a transient term and the curl of H as a source term that comes from the other coupled equation.

Probably I will have to adapt some code, because I want to implement a new (high-order) numerical scheme that I am doing for my master thesis. But I want to take advantage of the object-oriented structure of FiPy to make reusable and redistributable code.

Is there any example on how would this be done? I have looked at the documentation but I couldn't find something that I could translate properly on this.

[guyer]

I take it you mean 2nd order time discretization, not second time derivative $\partial^2 / \partial t^2$ ?

[wd15]

I take it you mean 2nd order time discretization, not second time derivative $\partial^2 / \partial t^2$ ?

No, I think that @yurigba is referring to second order derivative like in Maxwell's equations or the wave equation. I think you break that down into first order derivative with more equations to solve it though.

[wd15]

Overall, FiPy needs quite a lot of work to have a good FDTD method or to solve hyperbolic equations well. Not that I want to dissuade your interest. We'd really welcome the contribution, but I want to be realistic about the amount of work involved.

I started looking at Riemann solvers for FiPy. I got the first order version working, see this. Unfortunately, I never got to the point of merging this back into FiPy, which I now regret. It might be helpful to you though.

You might want to look into CLAWPACK, which deals only with hyperbolic problems. I used this to check that the Riemann branch was working correctly.

For FDTD method I played with MEEP a bit. My hope was to try and reimplement some of that in FiPy, but again work that was never finished.

I would certainly start off by solving the problem with a code that specializes in the equations of interest and then reverse engineer into FiPy if you still want to.

[guyer]

Ah, OK. But then you don't need curls, do you?

[yurigba]

That is correct @wd15. Thanks for the sources, they may be very helpful in the future.

I understand the picture better now. Due to the fact that finite volume method is not the preferred way for electromagnetic problems, it's understandable that there should be few options to do this kind of approach.

Ultimately, using a conservation law, the curl integrated over a surface will become a circulation integral due to Stokes' theorem. So in practice what I need is to approximate the circulation over a 2D volume based on vertex values. This is pretty much what I need to be able to calculate. Time leap of E fields in the cell center approximately equals circulation of H on vertex nodes.

As for the transient term, I don't need, from the beginning, to be able to discretize the second order time derivative. However, for a second order scheme the transient term looks mostly like a wave equation. So I wanted to apply this nice property somehow.

In either case, for what I need to do it will always be a lot of work, because I will have to build from the ground up the new method that I am working. I just need to know what tools that are available.

Even if there is not a way to do this, I suppose it should not be so hard to use the meshes and variables packages, right? Is there some sample code of a simple solver that uses these packages by updating manually the values of the class variable?

[guyer]

@yurigba FiPy does not store or solve for vertex values. It's a cell-centered finite volume code. Staggered grids are an option (and are how the FDTD codes I looked at ages ago worked), but there may be a cell-centered FVTD formulation that makes more sense. Abstractly, I think it should be possible to create Term objects that achieve the discretizations you want. If you get it working, these would be most welcome contributions.

@wd15 and I are both attending SciPy 2020 right now, so it'll be a few days before I, at least, have the bandwidth to walk you through our classes. fipy.meshes and fipy.variables were never written to work independently of the rest of FiPy, but it would interesting to try and we're certainly open to making adjustments to get them working for you.

Is there some sample code of a simple solver that uses these packages by updating manually the values of the class variable?

We provide dozens of examples that show different ways of manipulating Variable objects and building equations to .solve() them. I get the sense that's not what you're asking, though. Can you explain a bit more?

[wd15]

That is correct @wd15. Thanks for the sources, they may be very helpful in the future.

I understand the picture better now. Due to the fact that finite volume method is not the preferred way for electromagnetic problems, it's understandable that there should be few options to do this kind of approach.

I wouldn't say that necessarily. You can probably use FD, FV or FE to solve these. I think that FDTD type problems could easily be solved using the FV and the distinction is not that important.

Ultimately, using a conservation law, the curl integrated over a surface will become a circulation integral due to Stokes' theorem. So in practice what I need is to approximate the circulation over a 2D volume based on vertex values. This is pretty much what I need to be able to calculate. Time leap of E fields in the cell center approximately equals circulation of H on vertex nodes.

As @guyer said this may be possible with a cell centered approach.

Even if there is not a way to do this, I suppose it should not be so hard to use the meshes and variables packages, right? Is there some sample code of a simple solver that uses these packages by updating manually the values of the class variable?

Any variable object in FiPy can be updated with the standard assignment operator. It's better to keep a clear distinction between variables that depend on other variables and those that are updated manually. The mesh module most likely has the arrays required for the calculation. I'm more worried about higher order accuracy and resolving waves correctly which are required for these types of problems. Anyway, if you need more help we will try our best.

[yurigba]

@guyer, @wd15, your help is much appreciated. From my own experience developing numerical methods i guess that there is no good solver on hyperbolic coupled problems without staggered grids. And my master thesis goal to give some intuition on this.

There is a geometrical link between FDTD and FVTD methods in electromagnetics based on geometric algebra and there is already an article on that.

My objective is to build on this idea to build a better and more general unstructured solver, for engineering purposes.

I have tried to build it in MATLAB, but I had no knowledge about mesh generation and low skill on unstructured solvers, and I started looking for alternatives. I found that FiPy would suit this need.

To achieve this, I have thought to use FiPy like this:

from fipy import CellVariable
from fipy.meshes.uniformGrid2D import UniformGrid2D
from fipy.terms.sourceTerm import SourceTerm

nx = 20
ny = nx
dx = 1.
dy = dx
L = dx*nx

meshH = UniformGrid2D(nx = nx, ny = ny, dx = dx, dy = dy, origin = (0)) (the primary mesh)
meshE = UniformGrid2D(nx = nx-1, ny = ny-1, dx = dx, dy = dy, origin = (0.5)) (the dual mesh)

Ez = CellVariable(
    name = "electric field z component",
    mesh = meshE,
    value = 0., # Initial conditions array here
    rank = 0, # the z component of electric field, a scalar
)

H = CellVariable(
    name = "magnetic field",
    mesh = meshH,
    value = 0., # initial conditions array here
    rank = 1, # the x and y magnetic field components
)

(pseudocode: build two arrays VertexToCenter that links _meshE_ 
 cell centers to _meshH_ cell vertex and vice-versa)

for t in range(0, time steps):

    (pseudocode: look for vertex nodes in _meshE_,
     look for correspondent cell centers in VertexToCenter,
     update the E fields in CellVariable in _meshE_ based on 
     the correspondent cell center nodes in _meshH_ as a SourceTerm)
    (pseudocode: look for vertex nodes in _meshH_,
     look for correspondent cell centers in VertexToCenter,
     update the H fields in CellVariable in _meshH_ based on 
     the correspondent cell center nodes in _meshE_ as a SourceTerm)

    plot solution

The idea is that if I have access to the field variables, I can just build a function that returns the weights of each node in the mesh to the update equation and update the CellVariable manually.

@wd15 I will try to do the assignment, if this works it may be easy to do the rest.

[yurigba]

I have done these functions as I stated, and I could use the mesh and variables packages as standalone and it worked also with the viewers package. So that's something. However, I missed some features that could be of help in this.

There should be some equivalence in the sense of module implementation of the UniformGrid2D and other kind of meshes, in the sense that I would use the same handles for the meshes. I am referencing issue #706, which is exactly something that would be handy in further generalization of code. I didn't figure yet how to work with boundary nodes in a general setting (without using the particularity of a UniformGrid).

Also, something that would be useful and I could contribute if I understand properly the meshes package, is to manually assign the CellCenters of a mesh to the VertexNodes of another mesh, and adding the connectivity based on adjacent cells, thus creating a dual mesh generator. That would not only help the solution of electromagnetic problems but also would be better for navier-stokes equations solvers, which are feasible to do by using the same tools that I did.

Apart from that, I have successfully done a simulation of my scheme in cartesian grids by following @wd15 advice, and I look forward on using these packages on more complex problems (material interfaces, more complex boundaries), which I will soon be trying.

[guyer]

@yurigba the XxxGridND classes are certainly intended to be completely compatible with other Mesh classes, but there are unquestionably gaps in the API. I don't think #706 should be hard to fix. Let us know if you encounter other similar discrepancies.

I think your dual mesh generator idea should be feasible. A general Mesh is instantiated with vertexCoords, faceVertexIDs, and cellFaceIDs.

vertexCoords come from .cellCenters of your existing mesh.

Connectivity to determine faceVertexIDs and cellFaceIDs will come from .cellToCellIDs (I think).

[yurigba]

Thanks @guyer. I liked a lot the way that Mesh was built, it looks very self-contained.

I have tried the following (just for understanding of what to expect):

meshH = UniformGrid2D(nx = nx-1, ny = ny-1, dx = dx, dy = dy, origin = (0))

vc = meshH.vertexCoords
fvid = meshH.faceVertexIDs
cfid = meshH.cellFaceIDs

test = Mesh(vertexCoords=vc, faceVertexIDs=fvid, cellFaceIDs=cfid)

and got the following error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-30-7ddc14cfbb14> in <module>
      5 cfid = meshH.cellFaceIDs
      6 
----> 7 test = Mesh(vertexCoords=vc, faceVertexIDs=fvid, cellFaceIDs=cfid)

~/Documents/Projects/git/clifipy/fipy/meshes/mesh.py in __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, communicator, _RepresentationClass, _TopologyClass)
     53 
     54         self._setTopology()
---> 55         self._setGeometry(scaleLength = 1.)
     56 
     57     """

~/Documents/Projects/git/clifipy/fipy/meshes/mesh.py in _setGeometry(self, scaleLength)
    129         (self._internalCellDistances,
    130          self._cellDistanceVectors) = self._calcCellDistAndVec()
--> 131         self.faceNormals = self._calcFaceNormals()
    132         self._orientedFaceNormals = self._calcOrientedFaceNormals()
    133         self._cellVolumes = self._calcCellVolumes()

~/Documents/Projects/git/clifipy/fipy/meshes/mesh.py in _calcFaceNormals(self)
    195         faceVertexCoords = numerix.take(self.vertexCoords, faceVertexIDs, axis=1)
    196         t1 = faceVertexCoords[:, 1,:] - faceVertexCoords[:, 0,:]
--> 197         t2 = faceVertexCoords[:, 2,:] - faceVertexCoords[:, 1,:]
    198         norm = numerix.cross(t1, t2, axis=0)
    199         ## reordering norm's internal memory for inlining

IndexError: index 2 is out of bounds for axis 1 with size 2

I suspect this is another incompatibility of the classes. I'm doing these tests to better understand how faceVertexIDs and cellFaceIDs are built. Am I doing something wrong?

[wd15]

It's not accounting for 2D faces only having 2 attached vertices. That's not very good (edit: actually probably not intended to work like that). I'm not sure that calling Mesh directly is ever tested for anything, but 3D meshes or if it ever was intended to work for 2D. Both of the following workarounds are possible.

from fipy.meshes.uniformGrid3D import UniformGrid3D
from fipy.meshes.mesh import Mesh

meshH = UniformGrid3D(nx=10, ny=10, nz=1, dx=1.0, dy=1.0, dz=1.0)

vc = meshH.vertexCoords
fvid = meshH.faceVertexIDs
cfid = meshH.cellFaceIDs

test = Mesh(vertexCoords=vc, faceVertexIDs=fvid, cellFaceIDs=cfid)

or

from fipy.meshes.uniformGrid2D import UniformGrid2D
from fipy.meshes.mesh2D import Mesh2D

meshH = UniformGrid2D(nx=10, ny=10, dx=1.0, dy=1.0)

vc = meshH.vertexCoords
fvid = meshH.faceVertexIDs
cfid = meshH.cellFaceIDs

test = Mesh2D(vertexCoords=vc, faceVertexIDs=fvid, cellFaceIDs=cfid)