2D Advection-diffusion with Taylor-Green Vortex

import sys
import feelpp
import feelpp.toolboxes.core as tb
from feelpp.toolboxes.cfpdes import *

sys.argv = ["feelpp_cfpdes_Taylor_Green"]
e = feelpp.Environment(sys.argv,
                        opts=tb.toolboxes_options("coefficient-form-pdes", "cfpdes"),
                        config=feelpp.globalRepository("cfpdes-Taylor_Green"))

def generateGeometry(filename, dim=2, hsize=0.1):
    """create gmsh mesh

    Args:
        filename (str): name of the file
        dim (int): dimension of the mesh
        hsize (float): mesh size
    """
    geo = """SetFactory("OpenCASCADE");
    h = {};
    dim = {};
    """.format(hsize, dim)
    if dim == 2:
        geo += """
        Rectangle(1) = {-1, -1, 0, 2, 2, 0};
        Characteristic Length{ PointsOf{ Surface{1}; } } = h;
        Physical Curve("Gamma_D") = {1};
        Physical Curve("Gamma_N") = {2,3,4};
        Physical Surface("Omega") = {1};
        """
    elif dim == 3:
        geo += """
        Box(1) = {-1, -1, -1, 2, 2, 2};
        Characteristic Length{ PointsOf{ Volume{1}; } } = h;
        Physical Surface("Gamma_D") = {5};
        Physical Surface("Gamma_N") = {1,2,3,4,6};
        Physical Volume("Omega") = {1};
        """
    with open(filename, 'w') as f:
        # Write the string to the file
        f.write(geo)

def getMesh(filename, hsize=0.05, dim=2, verbose=False):
    """create mesh

    Args:
        filename (str): name of the file
        hsize (float): mesh size
        dim (int): dimension of the mesh
        verbose (bool): verbose mode
    """
    import os
    for ext in [".msh",".geo"]:
        f = os.path.splitext(filename)[0] + ext
        if os.path.exists(f):
            os.remove(f)
    if verbose:
        print(f"generate mesh {filename} with hsize={hsize} and dimension={dim}")
    generateGeometry(filename=filename, dim=dim, hsize=hsize)
    mesh = feelpp.load(feelpp.mesh(dim=dim, realdim=dim), filename, hsize)
    return mesh

Taylor_Green_json = lambda order,dim=2,mu=(1/30),name="c":  {
    "Name": "Taylor_Green",
    "ShortName": "Taylor_Green",
    "Models":
    {
        f"cfpdes-{dim}d":
        {
            "equations":"AdvDiff"
        },
        "AdvDiff":{
            "setup":{
                "unknown":{
                    "basis":f"Pch{order}",
                    "name":f"{name}",
                    "symbol":"c"
                },
                "coefficients":{
                    "d":"1.",
                    "c":f"{mu}",
                    "beta":"{sin(pi*x)*cos(pi*y),-cos(pi*x)*sin(pi*y)}:x:y" if dim==2 else "{sin(pi*x)*cos(pi*y)*cos(pi*z),-0.5*cos(pi*x)*sin(pi*y)*cos(pi*z),-0.5*cos(pi*x)*cos(pi*y)*sin(pi*z)}:x:y:z",
                    "f":"0."
                }
            }
        }
    },

    "Materials":
    {
        "Omega":
        {
            "markers":["Omega"]
        }
    },
    "BoundaryConditions":
    {
        "AdvDiff":
        {
            "Dirichlet":
            {
                "g":
                {
                    "markers":["Gamma_D"],
                    "expr":"0"
                }
            },
            "Neumann":{
                "m":
                {
                    "markers":["Gamma_N"],
                    "expr":"0"
                }
            }
        }
    },
    "InitialConditions":
    {
        "AdvDiff":
        {
            f"{name}":
            {
                "Expression":
                {
                    "wendland":
                    {
                        "markers":["Omega"],
                        "expr":"((1-2.5*sqrt((x+0.6)^2+(y+0.6)^2))>0)*(1-2.5*sqrt((x+0.6)^2+(y+0.6)^2))^3*(3*(2.5*sqrt((x+0.6)^2+(y+0.6)^2))+1)+((1-2.5*sqrt((x)^2+(y)^2))>0)*(1-2.5*sqrt((x)^2+(y)^2))^3*(3*(2.5*sqrt((x)^2+(y)^2))+1)+((1-2.5*sqrt((x-0.6)^2+(y-0.6)^2))>0)*(1-2.5*sqrt((x-0.6)^2+(y-0.6)^2))^3*(3*(2.5*sqrt((x-0.6)^2+(y-0.6)^2))+1):x:y" if dim==2 else "((1-2.5*sqrt((x+0.6)^2+(y+0.6)^2+(z+0.6)^2))>0)*(1-2.5*sqrt((x+0.6)^2+(y+0.6)^2+(z+0.6)^2))^3*(3*(2.5*sqrt((x+0.6)^2+(y+0.6)^2+(z+0.6)^2))+1)+((1-2.5*sqrt((x)^2+(y)^2+(z)^2))>0)*(1-2.5*sqrt((x)^2+(y)^2+(z)^2))^3*(3*(2.5*sqrt((x)^2+(y)^2+(z)^2))+1)+((1-2.5*sqrt((x-0.6)^2+(y-0.6)^2+(z-0.6)^2))>0)*(1-2.5*sqrt((x-0.6)^2+(y-0.6)^2+(z-0.6)^2))^3*(3*(2.5*sqrt((x-0.6)^2+(y-0.6)^2+(z-0.6)^2))+1):x:y:z"
                    }
                }
            }
        }
    },
    "PostProcess":
    {
        f"cfpdes-{dim}d":
        {
            "Exports":
            {
                "fields":["all"]
            }
        }
    }
}

Taylor_Green_cfg = lambda time_order=2, t0=0, T=0.2, dt=0.1: f"""
case.dimensions=2
[cfpdes]
#filename=$cfgdir/Taylor_Green.json
pc-type=gamg
ksp-converged-reason=
#verbose=1
#solver=Newton
snes-monitor=1
reuse-prec=1
#use-cst-matrix=0
#use-cst-vector=0
#snes-line-search-type=l2#basic
[cfpdes.AdvDiff]
#time-stepping=Theta
#stabilization=1
#stabilization.type=unusual-gls #supg#unusual-gls #gls
[cfpdes.AdvDiff.bdf]
order={time_order}

[ts]
time-initial={t0}
time-step={dt}
time-final={T}
restart.at-last-save=true
"""

SolveTG(0.04, Taylor_Green_json(order=2, dim=2, mu=1/30, name="c"), "Taylor_Green.cfg")