2D cylindrical coordinates

Hi all, hi Juan,

I have a question regarding projects based on 2D cylindrical symmetry.
I tried to convert the 2D capacitor example into a project with Cylindrical symmetry using the commands described in the manual.

I thought it would be sufficient to add

ds.set_parameter(name="raxis_zero", value=0) ds.set_parameter(name="raxis_variable", value="x")

ds.cylindrical_node_volume(device=device, region=region)

ds.cylindrical_edge_couple(device=device, region=region)

ds.cylindrical_surface_area(device=device, region=region)

as described in the manual in section 4.6, and then to change the NodeVolumes, EdgeCoupleVolumes and EdgeNodeVolumes with set_parameter() for all the used models.

The simulation runs without error, and there are messages indicating a change to a cylindrical coordinate system, but the result (e.g. the contact charge values) are still the same as if I stick to the usual cartesian coordinate system.

I would expect a change in these charge, since the area of the capacitor plates should change in the cylindrical case? Or should they be identical, and I am missing something in the geometry?

Can anyone give me a hint, what I might be doing wrong?

Here is the full code, based on the 2D capacitor example:

import devsim as ds

device="MyDevice"
region="MyRegion"

xmin=-0.5025
x1  =-0.5
x2  =-0.5
x3  =0.5
x4  =0.5
xmax=0.5025

ymin=0.0
y1  =0.1
y2  =0.2
y3  =1.2
y4  =1.3
ymax=5

ds.create_2d_mesh(mesh=device)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=ymin, ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=y1  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=y2  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=y3  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=y4  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="y", pos=ymax, ps=0.5)

device=device
region="air"

ds.add_2d_mesh_line(mesh=device, dir="x", pos=xmin, ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="x", pos=x1  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="x", pos=x2  , ps=0.05)
ds.add_2d_mesh_line(mesh=device, dir="x", pos=x3  , ps=0.05)
ds.add_2d_mesh_line(mesh=device, dir="x", pos=x4  , ps=0.1)
ds.add_2d_mesh_line(mesh=device, dir="x", pos=xmax, ps=0.1)

ds.add_2d_region(mesh=device, material="gas"  , region="air", yl=ymin, yh=ymax, xl=xmin, xh=xmax)
ds.add_2d_region(mesh=device, material="metal", region="m1" , yl=y1  , yh=y2  , xl=x1  , xh=x4)
ds.add_2d_region(mesh=device, material="metal", region="m2" , yl=y3  , yh=y4  , xl=x2  , xh=x3)

# must be air since contacts don't have any equations
ds.add_2d_contact(mesh=device, name="bot", region="air", yl=y1, yh=y2, xl=x1, xh=x4, material="metal")
ds.add_2d_contact(mesh=device, name="top", region="air", yl=y3, yh=y4, xl=x2, xh=x3, material="metal")
ds.finalize_mesh(mesh=device)
ds.create_device(mesh=device, device=device)



###
### Set parameters on the region
###
ds.set_parameter(device=device, region=region, name="Permittivity", value=3.9*8.85e-14)

###
### Create the Potential solution variable
###
ds.node_solution(device=device, region=region, name="Potential")
ds.set_parameter(device=device, region=region, 
                     name="Potential", value="CylindricalNodeVolume")
###
### Creates the Potential@n0 and Potential@n1 edge model
###
ds.edge_from_node_model(device=device, region=region, node_model="Potential")

###
### Electric field on each edge, as well as its derivatives with respect to
### the potential at each node
###
ds.edge_model(device=device, region=region, name="ElectricField",
           equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")
ds.set_parameter(device=device, region=region, 
                     name="ElectricField", value="CylindricalEdgeCouple")

ds.edge_model(device=device, region=region, name="ElectricField:Potential@n0",
           equation="EdgeInverseLength")
ds.set_parameter(device=device, region=region,
                     name="ElectricField:Potential@n0", value="CylindricalEdgeNodeVolume@n0")

ds.edge_model(device=device, region=region, name="ElectricField:Potential@n1",
           equation="-EdgeInverseLength")
ds.set_parameter(device=device, region=region, 
                     name="ElectricField:Potential@n1", value="CylindricalEdgeNodeVolume@n1")
###
### Model the D Field
###
ds.edge_model(device=device, region=region, name="DField",
           equation="Permittivity*ElectricField")
ds.set_parameter(device=device, region=region, 
                     name="DField", value="CylindricalEdgeCouple")

ds.edge_model(device=device, region=region, name="DField:Potential@n0",
           equation="diff(Permittivity*ElectricField, Potential@n0)")
ds.set_parameter(device=device, region=region,
                     name="DField:Potential@n0", value="CylindricalEdgeNodeVolume@n0")


ds.edge_model(device=device, region=region, name="DField:Potential@n1",
           equation="-DField:Potential@n0")
ds.set_parameter(device=device, region=region,
                     name="DField:Potential@n1", value="CylindricalEdgeNodeVolume@n0")

###
### Create the bulk equation
###
ds.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
         edge_model="DField", variable_update="default")


###
### Contact models and equations
###
for c in ("top", "bot"):
    ds.contact_node_model(device=device, contact=c, name="%s_bc" % c,
                       equation="Potential - %s_bias" % c)

    ds.contact_node_model(device=device, contact=c, name="%s_bc:Potential" % c,
                       equation="1")

    ds.contact_equation(device=device, contact=c, name="PotentialEquation",
                     variable_name="Potential",
                     node_model="%s_bc" % c, edge_charge_model="DField")

###
### Set the contact 
###
ds.set_parameter(device=device, name="top_bias", value=1.0e-0)
ds.set_parameter(device=device, name="bot_bias", value=0.0)



ds.edge_model(device=device, region="m1", name="ElectricField", equation="0")


ds.edge_model(device=device, region="m2", name="ElectricField", equation="0")


ds.node_model(device=device, region="m1", name="Potential", equation="bot_bias;")
ds.node_model(device=device, region="m2", name="Potential", equation="top_bias;")

ds.set_parameter(name="raxis_zero", value=0)
ds.set_parameter(name="raxis_variable", value="x")

ds.cylindrical_node_volume(device=device, region=region)
ds.cylindrical_edge_couple(device=device, region=region)
ds.cylindrical_surface_area(device=device, region=region)

#solve -type dc -absolute_error 1.0 -relative_error 1e-10 -maximum_iterations 100 -solver_type iterative
ds.solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=30, solver_type="direct")

ds.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
#%%
Q_top = ds.get_contact_charge(device=device, contact="top", equation="PotentialEquation")
Q_bot = ds.get_contact_charge(device=device, contact="bot", equation="PotentialEquation")
print('Top contact charge is: '+str(Q_top))
print('Bottom contact charge is: '+str(Q_bot))
#%%
ds.write_devices(file="cap2d.msh", type="devsim")
ds.write_devices(file="cap2d.dat", type="tecplot")
ds.write_devices(file="cap2d", type="vtk")

As always, thank you very much for any help.

Please see this example:

examples/bioapp1/bioapp1_2d.py

and related blog post:
Biosensor Applications in DEVSIM | DEVSIM

Specifically in:

examples/bioapp1/bioapp1_common.py

you will see this code after creating the models:

set_parameter(name="node_volume_model",          value="CylindricalNodeVolume")
set_parameter(name="edge_couple_model",          value="CylindricalEdgeCouple")
set_parameter(name="element_edge_couple_model",  value="ElementCylindricalEdgeCouple")
set_parameter(name="element_node0_volume_model", value="ElementCylindricalNodeVolume@en0")
set_parameter(name="element_node1_volume_model", value="ElementCylindricalNodeVolume@en1")

This description here may be helpful:
4. Equation and Models — DEVSIM Manual 2.0.0 documentation

Please let me know if you get your expected result. This feature has not been widely used.

Great, thank you for the additional references to the files, I will try it and come back to you later.

Hi @Juan,

the 2D coordinate system works well, I have just added the above functions to the 2D capacitor file like you said and like it is done in your referenced example. Sorry for the earlier confusion, I assumed wrongly that I have to change the name to the model-names in my file.

ds.set_parameter(device=device, name="raxis_variable", value="x")
ds.set_parameter(device=device, name="raxis_zero",     value=0)
ds.cylindrical_node_volume(device=device, region=region)
ds.cylindrical_edge_couple(device=device, region=region)

ds.set_parameter(name="node_volume_model",          value="CylindricalNodeVolume")
ds.set_parameter(name="edge_couple_model",          value="CylindricalEdgeCouple")
ds.set_parameter(name="element_edge_couple_model",  value="ElementCylindricalEdgeCouple")
ds.set_parameter(name="element_node0_volume_model", value="ElementCylindricalNodeVolume@en0")
ds.set_parameter(name="element_node1_volume_model", value="ElementCylindricalNodeVolume@en1")

Here are some results to confirm:
This is the structure layout of the plate capacitor, which I simulated in standard 2D cartesian coordinates and with cylindrical 2D coordinates with rotation axis at x = 0. The extension of the plates in x-direction is dx = 1 cm in both cases, the distance between the plate is 1 cm and the capacitor is filled with SiO2 with permittivity of 3.9. The applied bias is 1 V.

The different results in accumulated charge Q [C] in capacitor for different plate extension dx are displayed here for both cases, where the lines correspond to the expected Q from epsilon * A/d.
Since A_cyl = pi*dx^2 for the circular plates in the cylindrical case and A_cart = dx * 1 for the cartesian simulation, the increase is linear and quadratic, respectively.

cap2D_cartesian_vs_cylindrical

Just as a note: for the cylindrical case, the rotation axis should be at the side of a 2D cross-section - I thought if I would set it to x = 0 for the following case with symmetric dx of 1 cm in both directions, the resulting rotation would still be a circle with 1 cm radius.
The calculated Q in this case is however two times the Q of a 1 cm radius cylindrical capacitor in this orientation. Just wanted to note that if someone stumbles over it.

In conclusion, with accessing and visualizing the data directly in Python, as well as now understanding the 2D cylindrical symmetry correctly, I can now finally move on to the structures I actually want to simulate.

In case something interesting comes up, I will post it here.

Thank you very much for your patience and help, @Juan!

Hi @SinglePhotonGuy , thanks for your very interesting work. It is nice to see it works for your example. You are the first person besides me to use this feature.