diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 339ac716..04565a03 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -124,7 +124,7 @@ jobs: ${{ runner.os }}-test- ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@v1 - - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl poisson_transient.jl TopOptEMFocus.jl + - run: cd test/; julia --project=.. --color=yes --check-bounds=yes runtests.jl isotropic_damage.jl validation_DrWatson.jl interpolation_fe.jl transient_linear.jl transient_nonlinear.jl TopOptEMFocus.jl tutorials_mpi: name: TutorialsMPI ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} runs-on: ${{ matrix.os }} diff --git a/Project.toml b/Project.toml index 23ec8a46..3160d4e4 100644 --- a/Project.toml +++ b/Project.toml @@ -33,8 +33,8 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] -Gridap = "=0.17.19" -GridapDistributed = "=0.3" +Gridap = "0.18" +GridapDistributed = "0.4" GridapGmsh = "=0.7" GridapP4est = "=0.3" GridapPETSc = "=0.5" diff --git a/assets/poisson_transient/poisson_transient.gif b/assets/poisson_transient/poisson_transient.gif deleted file mode 100644 index 25dee0a9..00000000 Binary files a/assets/poisson_transient/poisson_transient.gif and /dev/null differ diff --git a/assets/transient_linear/result.gif b/assets/transient_linear/result.gif new file mode 100644 index 00000000..3772fd56 Binary files /dev/null and b/assets/transient_linear/result.gif differ diff --git a/assets/transient_nonlinear/result.gif b/assets/transient_nonlinear/result.gif new file mode 100644 index 00000000..c2b20a11 Binary files /dev/null and b/assets/transient_nonlinear/result.gif differ diff --git a/deps/build.jl b/deps/build.jl index eb90b6fb..9cb4104b 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -21,7 +21,8 @@ files = [ "On using DrWatson.jl"=>"validation_DrWatson.jl", "Interpolation of CellFields"=>"interpolation_fe.jl", "Poisson equation on parallel distributed-memory computers"=>"poisson_distributed.jl", - "Transient Poisson equation"=>"poisson_transient.jl", + "Transient Poisson equation"=>"transient_linear.jl", + "Transient nonlinear equation"=>"transient_nonlinear.jl", "Topology optimization"=>"TopOptEMFocus.jl", "Unfitted Poisson"=>"unfitted_poisson.jl"] diff --git a/docs/make.jl b/docs/make.jl index 9efaeb4f..cf642a0c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,7 +32,8 @@ nbviwer_logo = "https://img.shields.io/badge/show-nbviewer-579ACA.svg" for (i,(title,filename)) in enumerate(Tutorials.files) # Generate strings tutorial_prefix = string("t",@sprintf "%03d_" i) - tutorial_title = string("# # Tutorial ", i, ": ", title) + tutorial_id = replace(filename, " " => "_") + tutorial_title = string("# # [Tutorial ", i, ": ", title, "](@id ", tutorial_id, ")") tutorial_file = string(tutorial_prefix,splitext(filename)[1]) notebook_filename = string(tutorial_file, ".ipynb") binder_url = joinpath("@__BINDER_ROOT_URL__","notebooks", notebook_filename) diff --git a/src/TopOptEMFocus.jl b/src/TopOptEMFocus.jl index 9933e0de..2686996a 100644 --- a/src/TopOptEMFocus.jl +++ b/src/TopOptEMFocus.jl @@ -1,19 +1,19 @@ # In this tutorial, we will learn: -# +# # * How to apply the adjoint method for sensitivity analysis in Gridap # * How to do topology optimization in Gridap -# -# We recommend that you first read the [Electromagnetic scattering tutorial](https://gridap.github.io/Tutorials/dev/pages/t012_emscatter/#Tutorial-12:-Electromagnetic-scattering-in-2D-1) to make sure you understand the following points: -# +# +# We recommend that you first read the [Electromagnetic scattering tutorial](@ref emscatter.jl) to make sure you understand the following points: +# # * How to formulate the weak form for a 2d time-harmonic electromagnetic problem (a scalar Helmholtz equation) # * How to implement a perfectly matched layer (PML) to absorb outgoing waves # * How to impose periodic boundary conditions in Gridap # * How to discretize PDEs with complex-valued solutions -# +# # ## Problem statement -# +# # Consider the following optimization problem adapted from [Christiansen et al. (2020)](http://doi.org/10.1364/OE.28.004444): # We want to design a metallic (silver) nanoparticle to focus an incident $H_z$-polarized planewave on # a single spot, maximizing the electric-field intensity at this focal spot. The metallic @@ -22,95 +22,95 @@ # between a minimum radius $r_s = 10$nm (the minimum distance from the focal spot) and an outer # radius $r_d=100$nm. The computational cell is of height $H$ and length $L$, and we employ a perfectly matched layer (PML) thickness of $d_{pml}$ to implement outgoing (radiation) boundary conditions for this finite domain. # ![](../assets/TopOptEMFocus/Illustration.png) -# +# # The goal is find the arrangement of the silver material in the gray region that maximizes the |electric field|² at the center (the focal point). Every "pixel" in the gray region is effectively treated as a degree of freedom that can vary continuously between silver (shown in black below) and air (shown in white below). This is called density-based [topology optimization (TO)](https://en.wikipedia.org/wiki/Topology_optimization), and leads to a tractable optimization problem despite the huge number of parameters. A standard "projection" technique, described below, is used to "binarize" the structure by eventually forcing the material to be either silver or air almost everywhere. -# +# # ## Formulation -# -# From Maxwell's equations, considering a time-harmonic electromagnetic field polarized so that the electric field is in-plane and the magnetic field is out-of-plane (described by a scalar $H$ equal to the z-component), we can derive the governing equation of this problem in 2D (Helmholtz equation) [1]: -# +# +# From Maxwell's equations, considering a time-harmonic electromagnetic field polarized so that the electric field is in-plane and the magnetic field is out-of-plane (described by a scalar $H$ equal to the z-component), we can derive the governing equation of this problem in 2D (Helmholtz equation) [1]: +# # ```math -# \left[-\nabla\cdot\frac{1}{\varepsilon(x)}\nabla -k^2\mu(x)\right] H = f(x), +# \left[-\nabla\cdot\frac{1}{\varepsilon(x)}\nabla -k^2\mu(x)\right] H = f(x), # ``` -# +# # where $k=\omega/c$ is the wave number in free space and $f(x)$ is the source term (which corresponds to a magnetic current density in Maxwell's equations). -# +# # In order to simulate this scattering problem in a finite computational domain, we need outgoing (radiation) boundary conditions to prevent waves from reflecting back from the boundaries of the domain. We employ the well-known technique of "perfectly matched layers" (PML) [2], which are an an artificial absorbing layer adjacent to the boundaries that absorbs waves with minimal reflections (going to zero as the resolution increases). The "stretched-coordinate" formulation of PML correspond to a simple transformation of the PDE [3]: -# +# # ```math # \frac{\partial}{\partial x}\rightarrow \frac{1}{1+\mathrm{i}\sigma(u_x)/\omega}\frac{\partial}{\partial x}, # ``` -# +# # ```math # \frac{\partial}{\partial y}\rightarrow \frac{1}{1+\mathrm{i}\sigma(u_y)/\omega}\frac{\partial}{\partial y}, # ``` -# -# where $u_{x/y}$ is the depth into the PML, $\sigma$ is a profile function (here we chose $\sigma(u)=\sigma_0(u/d_{pml})^2$) and the $x$ and $y$ derivatives correspond PML layers at the $x$ and $y$ boundaries, respectively. Note that at a finite mesh resolution, PML reflects some waves, and the standard technique to mitigate this is to "turn on" the PML absorption gradually—in this case we use a quadratic profile. The amplitude $\sigma_0$ is chosen so that in the limit of infinite resolution the "round-trip" normal-incidence is some small number. -# -# Since PML absorbs all waves before they reach the boundary, the associated boundary condition can then be chosen arbitrarily. Here, the boundary conditions are Dirichlet (zero) on the top and bottom sides $\Gamma_D$ but periodic on the left ($\Gamma_L$) and right sides ($\Gamma_R$). The reason that we use periodic boundary conditions for the left and right side instead of Dirichlet boundary conditions is that we want to simulate a plane wave excitation, so we must choose boundary conditions that are satisfied by this incident wave. (Because of the anisotropic nature of PML, the PML layers at the $x$ boundaries do not disturb an incident planewave traveling purely in the $y$ direction.) -# -# Let $\mu(x)=1$ (materials at optical frequencies have negligible magnetic responses) and denote $\Lambda=\operatorname{diagm}(\Lambda_x,\Lambda_y)$ where $\Lambda_{x/y}=\frac{1}{1+\mathrm{i}\sigma(u_{x/y})/\omega}$. We can then formulate the problem as -# +# +# where $u_{x/y}$ is the depth into the PML, $\sigma$ is a profile function (here we chose $\sigma(u)=\sigma_0(u/d_{pml})^2$) and the $x$ and $y$ derivatives correspond PML layers at the $x$ and $y$ boundaries, respectively. Note that at a finite mesh resolution, PML reflects some waves, and the standard technique to mitigate this is to "turn on" the PML absorption gradually—in this case we use a quadratic profile. The amplitude $\sigma_0$ is chosen so that in the limit of infinite resolution the "round-trip" normal-incidence is some small number. +# +# Since PML absorbs all waves before they reach the boundary, the associated boundary condition can then be chosen arbitrarily. Here, the boundary conditions are Dirichlet (zero) on the top and bottom sides $\Gamma_D$ but periodic on the left ($\Gamma_L$) and right sides ($\Gamma_R$). The reason that we use periodic boundary conditions for the left and right side instead of Dirichlet boundary conditions is that we want to simulate a plane wave excitation, so we must choose boundary conditions that are satisfied by this incident wave. (Because of the anisotropic nature of PML, the PML layers at the $x$ boundaries do not disturb an incident planewave traveling purely in the $y$ direction.) +# +# Let $\mu(x)=1$ (materials at optical frequencies have negligible magnetic responses) and denote $\Lambda=\operatorname{diagm}(\Lambda_x,\Lambda_y)$ where $\Lambda_{x/y}=\frac{1}{1+\mathrm{i}\sigma(u_{x/y})/\omega}$. We can then formulate the problem as +# # ```math -# \left\{ \begin{aligned} +# \left\{ \begin{aligned} # \left[-\Lambda\nabla\cdot\frac{1}{\varepsilon(x)}\Lambda\nabla -k^2\right] H &= f(x) & \text{ in } \Omega,\\ # H&=0 & \text{ on } \Gamma_D,\\ # H|_{\Gamma_L}&=H|_{\Gamma_R},&\\ # \end{aligned}\right. # ``` -# +# # For convenience in the weak form and Julia implementation below, we represent $\Lambda$ as a vector given by the diagonal entries of the $2 \times 2$ scaling matrix, in which case $\Lambda\nabla$ becomes the elementwise product. -# +# # ## Topology Optimization -# +# # We use density-based topology optimization (TO) to maximize the electric field intensity at the center. In TO, every point in the design domain is a design degree of freedom that can vary continuously between air ($p=0$) and silver ($p=1$), which we discretize into a piece-wise constant parameter space $P$ for the design parameter $p\in [0,1]$. The material's electric permittivity ε is then given by: -# +# # ```math -# \varepsilon(p) = \left[n_{air}+p(n_{metal}-n_{air})\right]^2, +# \varepsilon(p) = \left[n_{air}+p(n_{metal}-n_{air})\right]^2, # ``` -# where $n_{air}=1$ and $n_{metal}$ are the refractive indices ($\sqrt{\varepsilon}$) of the air and metal, respectively. (It is tempting to simply linearly interpolate the permittivities ε, rather than the refractive indices, but this turns out to lead to artificial singularities in the case of metals where ε can pass through zero [4].) -# -# In practice, to avoid obtaining arbitrarily fine features as the spatial resolution is increased, one needs to regularize the problem with a minimum length-scale $r_f$ by generating a smoothed/filtered parameter function $p_f$. (Although this regularizes the problem, strictly speaking it does not impose a minimum feature size because of the nonlinear-projection step below. In practical applications, one imposes additional [manufacturing constraints](http://doi.org/10.1364/OE.431188) explicitly.) We perform the smoothing $p \to p_f$ by solving a simple "damped diffusion" PDE, also called a Helmholtz filter [5], for $p_f$ given the design variables $p$: +# where $n_{air}=1$ and $n_{metal}$ are the refractive indices ($\sqrt{\varepsilon}$) of the air and metal, respectively. (It is tempting to simply linearly interpolate the permittivities ε, rather than the refractive indices, but this turns out to lead to artificial singularities in the case of metals where ε can pass through zero [4].) +# +# In practice, to avoid obtaining arbitrarily fine features as the spatial resolution is increased, one needs to regularize the problem with a minimum length-scale $r_f$ by generating a smoothed/filtered parameter function $p_f$. (Although this regularizes the problem, strictly speaking it does not impose a minimum feature size because of the nonlinear-projection step below. In practical applications, one imposes additional [manufacturing constraints](http://doi.org/10.1364/OE.431188) explicitly.) We perform the smoothing $p \to p_f$ by solving a simple "damped diffusion" PDE, also called a Helmholtz filter [5], for $p_f$ given the design variables $p$: # ```math -# \begin{aligned} +# \begin{aligned} # -r_f^2\nabla^2p_f+p_f&=p\, ,\\ -# \left. \frac{\partial p_f}{\partial \vec{n}} \right\vert_{\partial\Omega_D} & =0 . -# \end{aligned} +# \left. \frac{\partial p_f}{\partial \vec{n}} \right\vert_{\partial\Omega_D} & =0 . +# \end{aligned} # ``` # -# We choose a filter radius $r_f=R_f/(2\sqrt{3})$ where $R_f=10$ nm, in order to match a published result (using a slightly different filtering scheme) for comparison [6]. -# -# Next, we apply a smoothed threshold projection to the intermediate variable $p_f$ to obtain a "binarized" density parameter $p_t$ that tends towards values of $0$ or $1$ almost everywhere [6] as the steepness $\beta$ of the thresholding is increased: +# We choose a filter radius $r_f=R_f/(2\sqrt{3})$ where $R_f=10$ nm, in order to match a published result (using a slightly different filtering scheme) for comparison [6]. +# +# Next, we apply a smoothed threshold projection to the intermediate variable $p_f$ to obtain a "binarized" density parameter $p_t$ that tends towards values of $0$ or $1$ almost everywhere [6] as the steepness $\beta$ of the thresholding is increased: # ```math -# p_t = \frac{\tanh(\beta\eta)+\tanh\left[\beta(p_f-\eta)\right]}{\tanh(\beta\eta)+\tanh\left[\beta(1-\eta)\right]}. +# p_t = \frac{\tanh(\beta\eta)+\tanh\left[\beta(p_f-\eta)\right]}{\tanh(\beta\eta)+\tanh\left[\beta(1-\eta)\right]}. # ``` -# Note that as $\beta\to\infty$, this threshold procedure goes to a step function, which would make the optimization problem non-differentiable. In consequence, the standard approach is to gradually increase $\beta$ to slowly binarize the design as the optimization progresses [6]. We will show how this is done below. -# +# Note that as $\beta\to\infty$, this threshold procedure goes to a step function, which would make the optimization problem non-differentiable. In consequence, the standard approach is to gradually increase $\beta$ to slowly binarize the design as the optimization progresses [6]. We will show how this is done below. +# # Finally, we replace $p$ with the filtered and thresholded $p_t$ in the ε interpolation formula from above: -# +# # ```math # \varepsilon(p_t) = \left[n_{air}+p_t(n_{metal}-n_{air})\right]^2, # ``` # This is the quantity that will be used for the $1/\varepsilon(x)$ coefficient in our Helmholtz PDE. -# +# # ## Weak form -# +# # Now we derive the weak form of the PML Helmholtz PDE above. After integration by parts (in which the boundary term vanishes), we obtain: -# +# # ```math # a(u,v,p) = \int_\Omega \left[\nabla(\Lambda v)\cdot\frac{1}{\varepsilon(p)}\Lambda\nabla u-k^2uv\right]\mathrm{d}\Omega, # ``` -# +# # ```math -# b(v) = \int_\Omega vf\mathrm{d}\Omega. +# b(v) = \int_\Omega vf\mathrm{d}\Omega. # ``` # Notice that the $\nabla(\Lambda v)$ is an element-wise "product" of two vectors $\nabla$ and $\Lambda v$. -# +# # ## Setup -# -# We import the packages that will be used, define the geometry and physics parameters. -# +# +# We import the packages that will be used, define the geometry and physics parameters. +# using Gridap, Gridap.Geometry, Gridap.Fields, GridapGmsh λ = 532 # Wavelength (nm) @@ -125,16 +125,16 @@ n_air = 1 # Air refractive index k = 2*π/λ # Wavenumber (nm^-1) # ## Discrete Model -# -# We import the model from the `RecCirGeometry.msh` mesh file using the `GmshDiscreteModel` function defined in `GridapGmsh`. The mesh file is created with GMSH in Julia (see the file ../assets/TopOptEMFocus/MeshGenerator.jl). Note that this mesh file already specifies periodic boundaries for the left and right sides, which will cause Gridap to implement periodic boundary conditions. Also, the center smallest-distance circle region is labeled with `Center` and the annular design region is labeled with `Design` in the mesh file. -# +# +# We import the model from the `RecCirGeometry.msh` mesh file using the `GmshDiscreteModel` function defined in `GridapGmsh`. The mesh file is created with GMSH in Julia (see the file ../assets/TopOptEMFocus/MeshGenerator.jl). Note that this mesh file already specifies periodic boundaries for the left and right sides, which will cause Gridap to implement periodic boundary conditions. Also, the center smallest-distance circle region is labeled with `Center` and the annular design region is labeled with `Design` in the mesh file. +# model = GmshDiscreteModel("../models/RecCirGeometry.msh") # ## FE spaces for the magnetic field -# +# # We use the first-order Lagrange finite-element basis functions. The Dirichlet edges are labeled as `DirichletEdges` in the mesh file. Since our problem involves complex numbers (because of the PML and the complex metal refractive index), we need to specify the `vector_type` as `Vector{ComplexF64}`. -# +# order = 1 reffe = ReferenceFE(lagrangian, Float64, order) @@ -142,19 +142,19 @@ V = TestFESpace(model, reffe, dirichlet_tags = ["DirichletEdges", "DirichletNode U = V # mathematically equivalent to TrialFESpace(V,0) # ## Numerical integration -# +# # We construct the triangulation and a second-order Gaussian quadrature scheme for assembling the finite-element matrix from the weak form. Note that we create a boundary triangulation from a `Source` tag for the line excitation, which is a convenient and accurate way to produce an incident planewave. (We could have alternatively devised a corresponding current source, e.g. using a finite-width delta-function approximation.) -# +# degree = 2 Ω = Triangulation(model) dΩ = Measure(Ω, degree) -Γ_s = BoundaryTriangulation(model; tags = ["Source"]) # Source line +Γ_s = BoundaryTriangulation(model; tags = ["Source"]) # Source line dΓ_s = Measure(Γ_s, degree) # We also want to construct quadrature meshes for the numerical integration over two subsets of the computational cell: the design domain (annulus) $\Omega_d$ and the central "hole" $\Omega_c$ surrounding the focal point. The former is used to localize the design optimization to $\Omega_d$, and the latter is used to define the objective function (which only depends on the field at the center). -# +# Ω_d = Triangulation(model, tags="Design") dΩ_d = Measure(Ω_d, degree) @@ -164,9 +164,9 @@ dΩ_c = Measure(Ω_c, degree) # ## FE spaces for the design parameters -# +# # As discussed above, we need a piece-wise constant design parameter space $P$ that is defined in the design domain, this is achieved by a zero-order lagrangian. The number of design parameters is then the number of cells in the design region. -# +# p_reffe = ReferenceFE(lagrangian, Float64, 0) Q = TestFESpace(Ω_d, p_reffe, vector_type = Vector{Float64}) @@ -174,25 +174,25 @@ P = Q np = num_free_dofs(P) # Number of cells in design region (number of design parameters) # Note that this over 70k design parameters, which is large but not huge by modern standards. To optimize so many design parameters, the key point is how to compute the gradients to those parameters efficiently. -# +# # Also, we need a first-order lagrangian function space $P_f$ for the filtered parameters $p_f$ since the zero-order lagrangian always produces zero derivatives. (Note that $p_f$ and $p_t$ share the same function space since the latter is only a projection of the previous one.) -# +# pf_reffe = ReferenceFE(lagrangian, Float64, 1) Qf = TestFESpace(Ω_d, pf_reffe, vector_type = Vector{Float64}) Pf = Qf -# Finally, we pack up every thing related to Gridap as a named tuple called `fem_params`. This is because we want to pass those as local parameters to the optimization functions later, instead of making them as global parameters. -# +# Finally, we pack up every thing related to Gridap as a named tuple called `fem_params`. This is because we want to pass those as local parameters to the optimization functions later, instead of making them as global parameters. +# fem_params = (; V, U, Q, P, Qf, Pf, np, Ω, dΩ, dΩ_d, dΩ_c, dΓ_s) # ## PML formulation -# -# First we pack up all physical parameters as a structure call `phys`. Then we define a `s_PML` function: $s(x)=1+\mathrm{i}\sigma(u)/\omega,$ and its derivative `ds_PML`. The parameter `LHp` and `LHn` indicates the size of the inner boundary of the PML regions. Finally, we create a function-like object `Λ` that returns the PML factors and define its derivative in Gridap. -# +# +# First we pack up all physical parameters as a structure call `phys`. Then we define a `s_PML` function: $s(x)=1+\mathrm{i}\sigma(u)/\omega,$ and its derivative `ds_PML`. The parameter `LHp` and `LHn` indicates the size of the inner boundary of the PML regions. Finally, we create a function-like object `Λ` that returns the PML factors and define its derivative in Gridap. +# # Note that here we are defining a "callable object" of type `Λ` that encapsulates all of the PML parameters. This is convenient, both because we can pass lots of parameters around easily and also because we can define additional methods on `Λ`, e.g. to express the `∇(Λv)` operation. -# +# R = 1e-10 LHp=(L/2, h1+h2) # Start of PML for x,y > 0 @@ -229,13 +229,13 @@ Fields.∇(Λf::Λ) = x -> TensorValue{2, 2, ComplexF64}(-(Λf(x)[1])^2 * ds_PML # ## Filter and threshold -# +# # Here we use the filter and threshold discussed above. The parameters for the filter and threshold are extracted from Ref [6]. Note that every integral in the filter is only defined on $\Omega_d$ -# +# r = 5/sqrt(3) # Filter radius β = 32.0 # β∈[1,∞], threshold sharpness -η = 0.5 # η∈[0,1], threshold center +η = 0.5 # η∈[0,1], threshold center a_f(r, u, v) = r^2 * (∇(v) ⋅ ∇(u)) @@ -249,14 +249,14 @@ function Filter(p0; r, fem_params) end function Threshold(pfh; β, η) - return ((tanh(β * η) + tanh(β * (pfh - η))) / (tanh(β * η) + tanh(β * (1.0 - η)))) + return ((tanh(β * η) + tanh(β * (pfh - η))) / (tanh(β * η) + tanh(β * (1.0 - η)))) end # ## Weak form -# -# We notice that the design parameters only affect the weak form in the design domain and the PML does not affect the design domain, we can then make things simpler by dividing the weak form to a base integral that contains the whole computation cell and an additional integral on the design domain. We also make a LU factorization on the final Maxwell operator matrix $A$ since it will only be used to solve for linear equations. -# +# +# We notice that the design parameters only affect the weak form in the design domain and the PML does not affect the design domain, we can then make things simpler by dividing the weak form to a base integral that contains the whole computation cell and an additional integral on the design domain. We also make a LU factorization on the final Maxwell operator matrix $A$ since it will only be used to solve for linear equations. +# using LinearAlgebra ξd(p, n_air, n_metal)= 1 / (n_air + (n_metal - n_air) * p)^2 - 1 / n_air^2 # in the design region @@ -273,69 +273,69 @@ function MatrixA(pth; phys_params, fem_params) end # ## Solve for plane wave incident -# +# # The plane wave source `b_vec` can be simply assembled by a uniform integral over the source line, and the magnetic field vector `u_vec` can then be solved by a simple linear equation -# +# p0 = zeros(fem_params.np) # Here we make p=0 everywhere just for illustration purpose pf_vec = Filter(p0;r, fem_params) pfh = FEFunction(fem_params.Pf, pf_vec) pth = (pf -> Threshold(pf; β, η)) ∘ pfh -A_mat = MatrixA(pth; phys_params, fem_params) +A_mat = MatrixA(pth; phys_params, fem_params) b_vec = assemble_vector(v->(∫(v)fem_params.dΓ_s), fem_params.V) u_vec = A_mat \ b_vec uh = FEFunction(fem_params.U, u_vec) # ## Objective -# -# The problem is maximizing the electric field intensity at the center. Recall that the electric field can be retrieved from the magnetic field by +# +# The problem is maximizing the electric field intensity at the center. Recall that the electric field can be retrieved from the magnetic field by # # ```math # \mathbf{E}(\mathbf{x})=\frac{\mathrm{i}}{\omega\varepsilon(\mathbf{x})}\nabla\times\mathbf{H}(\mathbf{x}), -# ``` -# and our objective is the field intensity at center $\vert\mathbf{E}(\mathbf{x}_0)\vert^2$. -# -# In the 2D formulation, this objective can be simplified to +# ``` +# and our objective is the field intensity at center $\vert\mathbf{E}(\mathbf{x}_0)\vert^2$. +# +# In the 2D formulation, this objective can be simplified to # ```math # g = \int \vert\nabla H\vert^2\delta(x-x_0)\mathrm{d}\Omega = u^\dagger Ou, # ``` -# where $u$ is the magnetic field vector and +# where $u$ is the magnetic field vector and # ```math # O = \int (\nabla \hat{v}\cdot\nabla\hat{u})\delta(x-x_0)\mathrm{d}\Omega, # ``` -# with $\hat{v}$ and $\hat{u}$ are the finite element basis functions. -# -# In practice, the delta function can be approximated by a concentrated Gaussian function. Note that we use `dΩ_c` here in order to reduce computation costs. -# +# with $\hat{v}$ and $\hat{u}$ are the finite element basis functions. +# +# In practice, the delta function can be approximated by a concentrated Gaussian function. Note that we use `dΩ_c` here in order to reduce computation costs. +# function MatrixOf(fem_params) x0 = VectorValue(0,300) # Position of the field to be optimized - δ = 1 + δ = 1 return assemble_matrix(fem_params.U, fem_params.V) do u, v ∫((x->(1/(2*π)*exp(-norm(x - x0)^2 / 2 / δ^2))) * (∇(u) ⋅ ∇(v)) )fem_params.dΩ_c end end # ## Optimization with adjoint method -# -# Now that we have our objective to optimize, the next step is to find out the derivative to the design parameter $p$ in order to apply a gradient-based optimization algorithm. We will be using `ChainRulesCore` and `Zygote` packages. -# +# +# Now that we have our objective to optimize, the next step is to find out the derivative to the design parameter $p$ in order to apply a gradient-based optimization algorithm. We will be using `ChainRulesCore` and `Zygote` packages. +# using ChainRulesCore, Zygote import ChainRulesCore: rrule NO_FIELDS = ZeroTangent() -# Recall that our objective is $g=u^\dagger Ou$ and only $u=A(p)^{-1} b$ depends on the design parameters. The derivative of $g$ with respect to $p_t$ can be obtained via [adjoint method](https://math.mit.edu/~stevenj/18.336/adjoint.pdf): +# Recall that our objective is $g=u^\dagger Ou$ and only $u=A(p)^{-1} b$ depends on the design parameters. The derivative of $g$ with respect to $p_t$ can be obtained via [adjoint method](https://math.mit.edu/~stevenj/18.336/adjoint.pdf): # ```math -# \frac{\mathrm{d} g}{\mathrm{d}p_t}= -2\Re\left[w^\dagger\left(\frac{\mathrm{d}A}{\mathrm{d}p_t}u\right)\right], -# ``` -# where $w$ comes from the adjoint solve $A^\dagger w = Ou$. The final derivative with respect to $p$ can then be obtained via chain rules: +# \frac{\mathrm{d} g}{\mathrm{d}p_t}= -2\Re\left[w^\dagger\left(\frac{\mathrm{d}A}{\mathrm{d}p_t}u\right)\right], +# ``` +# where $w$ comes from the adjoint solve $A^\dagger w = Ou$. The final derivative with respect to $p$ can then be obtained via chain rules: # ```math -# \frac{\mathrm{d} g}{\mathrm{d}p}=\frac{\mathrm{d} g}{\mathrm{d}p_t}\cdot\frac{\mathrm{d} p_t}{\mathrm{d}p_f}\cdot\frac{\mathrm{d} p_f}{\mathrm{d}p} +# \frac{\mathrm{d} g}{\mathrm{d}p}=\frac{\mathrm{d} g}{\mathrm{d}p_t}\cdot\frac{\mathrm{d} p_t}{\mathrm{d}p_f}\cdot\frac{\mathrm{d} p_f}{\mathrm{d}p} # ``` -# +# # First we define some relative derivative functions: -# +# Dptdpf(pf, β, η) = β * (1.0 - tanh(β * (pf - η))^2) / (tanh(β * η) + tanh(β * (1.0 - η))) @@ -345,7 +345,7 @@ DAdpf(u, v, pfh; phys_params, β, η) = ((p -> Dξdpf(p, phys_params.n_air, phys # Then we create a function `gf_pf` that depends directly on $p_f$ and write out the derivative using adjoint method formula. Note that the threshold chainrule is already implemented in the functions above. -# +# function gf_pf(pf_vec; β, η, phys_params, fem_params) pfh = FEFunction(fem_params.Pf, pf_vec) @@ -353,7 +353,7 @@ function gf_pf(pf_vec; β, η, phys_params, fem_params) A_mat = MatrixA(pth; phys_params, fem_params) b_vec = assemble_vector(v->(∫(v)fem_params.dΓ_s), fem_params.V) u_vec = A_mat \ b_vec - + O_mat = MatrixOf(fem_params) real(u_vec' * O_mat * u_vec) end @@ -372,18 +372,18 @@ function Dgfdpf(pf_vec; β, η, phys_params, fem_params) b_vec = assemble_vector(v->(∫(v)fem_params.dΓ_s), fem_params.V) u_vec = A_mat \ b_vec O_mat = MatrixOf(fem_params) - + uh = FEFunction(fem_params.U, u_vec) w_vec = A_mat' \ (O_mat * u_vec) - wconjh = FEFunction(fem_params.U, conj(w_vec)) - + wconjh = FEFunction(fem_params.U, conj(w_vec)) + l_temp(dp) = ∫(real(-2 * DAdpf(uh, wconjh, pfh; phys_params, β, η)) * dp)fem_params.dΩ_d dgfdpf = assemble_vector(l_temp, fem_params.Pf) return dgfdpf end # Next we define the relation between $p_f$ and $p$, and obtain the derivative of the filter by again applying an adjoint method: -# +# function pf_p0(p0; r, fem_params) pf_vec = Filter(p0; r, fem_params) @@ -408,7 +408,7 @@ function Dgdp(dgdpf; r, fem_params) end # Finally, we pack up into a single function that takes `p` and returns our objective function, and which can optionally take a `grad` vector into which the gradient (computed by Zygote by composing our rules above) can be written in-place (as required for use in the NLopt optimization package). We also optionally record the value of the objective function from every call in order to save a record of the optimization process. -# +# function gf_p(p0::Vector; r, β, η, phys_params, fem_params) pf_vec = pf_p0(p0; r, fem_params) @@ -428,7 +428,7 @@ function gf_p(p0::Vector, grad::Vector; r, β, η, phys_params, fem_params) end # Using the following codes, we can check if we can get the derivatives correctly from the adjoint method by comparing it with the finite difference results. -# +# p0 = rand(fem_params.np) δp = rand(fem_params.np)*1e-8 @@ -439,9 +439,9 @@ g1 = gf_p(p0+δp, []; r, β, η, phys_params, fem_params) g1-g0, grad'*δp # ## Optimization with NLopt -# +# # Now we use NLopt.jl package to implement the MMA algorithm for optimization. Note that we start with $\beta=8$ and then gradually increase it to $\beta=32$ in consistent with Ref. [6]. -# +# using NLopt @@ -473,9 +473,9 @@ end @show g_opt # ## Results and plot -# +# # We use the CairoMakie.jl and GridapMakie.jl packages to plot the field as well as the optimized shape. Note that there might be multiple local optima for this problem, so different initial guesses might result in different optimized shapes. -# +# using CairoMakie, GridapMakie p0 = p_opt @@ -498,13 +498,13 @@ save("shape.png", fig) # ![](../assets/TopOptEMFocus/shape.png) -# +# -# For the electric field, recall that $\vert E\vert^2\sim\vert \frac{1}{\epsilon}\nabla H\vert^2$, the factor 2 below comes from the amplitude compared to the incident plane wave. We can see that the optimized shapes are very similar to the optimized shape in Ref. [6], proving our results. -# +# For the electric field, recall that $\vert E\vert^2\sim\vert \frac{1}{\epsilon}\nabla H\vert^2$, the factor 2 below comes from the amplitude compared to the incident plane wave. We can see that the optimized shapes are very similar to the optimized shape in Ref. [6], proving our results. +# maxe = 30 # Maximum electric field magnitude compared to the incident plane wave -e1=abs2(phys_params.n_air^2) +e1=abs2(phys_params.n_air^2) e2=abs2(phys_params.n_metal^2) fig, ax, plt = plot(fem_params.Ω, 2*(sqrt∘(abs((conj(∇(uh)) ⋅ ∇(uh))/(CellField(e1,fem_params.Ω) + (e2 - e1) * pth)))), colormap = :hot, colorrange=(0, maxe)) @@ -515,20 +515,19 @@ limits!(ax, -rplot, rplot, (h1)/2-rplot, (h1)/2+rplot) save("Field.png", fig) # ![](../assets/TopOptEMFocus/Field.png) -# +# # ## References -# +# # [1] [Wikipedia: Electromagnetic wave equation](https://en.wikipedia.org/wiki/Electromagnetic_wave_equation) -# +# # [2] [Wikipedia: Perfectly matched layer](https://en.wikipedia.org/wiki/Perfectly_matched_layer) -# +# # [3] A. Oskooi and S. G. Johnson, “[Distinguishing correct from incorrect PML proposals and a corrected unsplit PML for anisotropic, dispersive media](http://math.mit.edu/~stevenj/papers/OskooiJo11.pdf),” Journal of Computational Physics, vol. 230, pp. 2369–2377, April 2011. -# +# # [4] R. E. Christiansen, J. Vester-Petersen, S.P. Madsen, and O. Sigmund, “[A non-linear material interpolation for design of metallic nano-particles using topology optimization](https://pure.au.dk/portal/files/206950673/1_s2.0_S0045782518304328_main_1_.pdf),” Computer Methods in Applied Mechanics and Engineering , vol. 343, pp. 23–39, January 2019. -# +# # [5] B. S. Lazarov and O. Sigmund, "[Filters in topology optimization based on Helmholtz-type differential equations](https://en.wikipedia.org/wiki/Jacobi%E2%80%93Anger_expansion)", International Journal for Numerical Methods in Engineering, vol. 86, pp. 765-781, December 2010. -# +# # [6] R.E. Christiansen, J. Michon, M. Benzaouia, O. Sigmund, and S.G. Johnson, "[Inverse design of nanoparticles for enhanced Raman scattering](https://opg.optica.org/oe/fulltext.cfm?uri=oe-28-4-4444&id=426514)," Optical Express, vol. 28, pp. 4444-4462, February 2020. -# - +# diff --git a/src/interpolation_fe.jl b/src/interpolation_fe.jl index e1d70768..36c72678 100644 --- a/src/interpolation_fe.jl +++ b/src/interpolation_fe.jl @@ -145,7 +145,7 @@ writevtk(get_triangulation(gₕ), "target", cellfields=["gₕ"=>gₕ]) # functions in Raviart-Thomas spaces are vector-valued. The degrees of # freedom of the RT spaces are fluxes of the function across the edge # of the element. Refer to the -# [tutorial](https://gridap.github.io/Tutorials/dev/pages/t007_darcy/) +# [tutorial](@ref darcy.jl) # on Darcy equation with RT for more information on the RT # elements. diff --git a/src/poisson_transient.jl b/src/poisson_transient.jl deleted file mode 100644 index 7a473d3c..00000000 --- a/src/poisson_transient.jl +++ /dev/null @@ -1,119 +0,0 @@ -# ## Introduction - -# In this tutorial we will learn how to use [`GridapODEs.jl`](https://github.com/gridap/GridapODEs.jl) for approximating transient PDEs by using time marching schemes (method of lines). We consider the *heat equation*, a.k.a. the transient Poisson equation. - -# We will focus on the time discretization on the equations, assuming that the reader is familiar with the `Gridap` API for spatial finite element discretizations. See, e.g., [tutorial 1](https://gridap.github.io/Tutorials/stable/pages/t001_poisson/) for more details. - -# ## Problem statement - -# We solve the heat equation in a 2-dimensional domain $\Omega$, the unit square, with Homogeneous Dirichlet boundaries on the whole boundary $\partial \Omega$. We consider a time-dependent conductivity $\kappa(t)=1.0 + 0.95\sin(2\pi t)$, a time-dependent volumetric forcing term $f(t) = \sin(\pi t)$ and a constant Homogeneous boundary condition $g = 0.0$. The initial solution is $u(x,0) = u_0 = 0$. With these definitions, the strong form of the problem reads: - -# ```math -# \left\lbrace -# \begin{aligned} -# \frac{\partial u(t)}{\partial t} -\kappa(t)\Delta u(t) = f(t) \ &\text{in} \ \Omega,\\ -# u(t) = 0 \ &\text{on}\ \Gamma_{\rm D},\\ -# u(0) = 0 \ &\text{in}\ \Omega\\ -# \end{aligned} -# \right. -# ``` - -# The weak form of the problem reads: find $u(t)\in U_g(t)$ such that - -# ```math -# m(t,u,v) + a(t,u,v) = b(t,v)\quad \forall v\in \ V -# ``` - -# Note that $U_g(t)$ is a transient FE space, in the sense that Dirichlet boundary value of functions in $U_g$ _can_ change in time (even though this is not the case in this tutorial). The definition of $m(u,v)$, $a(u,v)$ and $b(v)$ is as follows: - -# ```math -# \begin{aligned} -# m(t,u,v) = \int_\Omega v\frac{\partial u}{\partial t} d\Omega, \\ -# a(t,u,v) = \int_\Omega \kappa(t) \nabla v\cdot \nabla u d\Omega, \\ -# b(t,v) = \int_\Omega v\ f(t) d\Omega -# \end{aligned} -# ``` - -# ## Discrete model and Triangulation - -# As usual, let us first load `Gridap`. -using Gridap - -# First, we define the `DiscreteModel` and the `Triangulation`. More details on this can be found in [tutorial 2](https://gridap.github.io/Tutorials/stable/pages/t002_validation/). - - -𝒯 = CartesianDiscreteModel((0,1,0,1),(20,20)) -Ω = Interior(𝒯) -dΩ = Measure(Ω,2) - -# ## FE space - -# In this tutorial we will use linear Lagrangian Finite Elements. -refFE = ReferenceFE(lagrangian,Float64,1) - -# The space of test functions is constant in time and is defined in steady problems: -V = TestFESpace(𝒯,refFE,dirichlet_tags="boundary") - -# The trial space is now a `TransientTrialFESpace`, which is constructed from a `TestFESpace` and a function (or vector of functions) for the Dirichlet boundary condition/s. In that case, the boundary condition function is a time-independent constant, but it could also be a time-dependent field depending on the coordinates $x$ and time $t$. -g(x,t::Real) = 0.0 -g(t::Real) = x -> g(x,t) -U = TransientTrialFESpace(V,g) - -# ## Weak form - -# The weak form of the problem follows the same structure as other `Gridap` tutorials, where we define the bilinear and linear forms to define the `FEOperator`. In this case we need to deal with time-dependent quantities and with the presence of time derivatives. The former is handled by passing the time, $t$, as an additional argument to the form, i.e. $a(t,u,v)$. The latter is defined using the time derivative operator `∂t`. - -# The most general way of constructing a transient FE operator is by using the `TransientFEOperator` function, which receives a residual, a Jacobian with respect to the unknown and a Jacobian with respect to the time derivative. -κ(t) = 1.0 + 0.95*sin(2π*t) -f(t) = sin(π*t) -res(t,u,v) = ∫( ∂t(u)*v + κ(t)*(∇(u)⋅∇(v)) - f(t)*v )dΩ -jac(t,u,du,v) = ∫( κ(t)*(∇(du)⋅∇(v)) )dΩ -jac_t(t,u,duₜ,v) = ∫( duₜ*v )dΩ -op = TransientFEOperator(res,jac,jac_t,U,V) - -# We can also take advantage of automatic differentiation techniques to compute both Jacobians and use the `TransientFEOperator` function sending just the residual. -op_AD = TransientFEOperator(res,U,V) - -# Alternatively, we can exploit the fact that the problem is linear and use the transient Affine FE operator signature `TransientAffineFEOperator`. In that case, we send a form for the mass contribution, $m$, a form for the stiffness contribution, $a$, and the forcing term, $b$. -m(t,u,v) = ∫( u*v )dΩ -a(t,u,v) = ∫( κ(t)*(∇(u)⋅∇(v)) )dΩ -b(t,v) = ∫( f(t)*v )dΩ -op_Af = TransientAffineFEOperator(m,a,b,U,V) - -# ### Alternative FE operator definitions - -# For time-dependent problems with constant coefficients, which is not the case of this tutorial, one could use the optimized operator `TransientConstantMatrixFEOperator`, which assumes that the matrix contributions ($m$ and $a$) are time-independent. That is: -m₀(u,v) = ∫( u*v )dΩ -a₀(u,v) = ∫( κ(0.0)*(∇(u)⋅∇(v)) )dΩ -op_CM = TransientConstantMatrixFEOperator(m₀,a₀,b,U,V) - -# Going further, if we had a problem with constant forcing term, i.e. constant force and constant boundary conditions, we could have used the `TransientConstantFEOperator`. In that case the linear form is also time-independent. -b₀(v) = ∫( f(0.0)*v )dΩ -op_C = TransientConstantFEOperator(m₀,a₀,b₀,U,V) - -# ## Transient solver - -# Once we have the FE operator defined, we proceed with the definition of the transient solver. First, we define a linear solver to be used at each time step. Here we use the `LUSolver`, but other choices are possible. -linear_solver = LUSolver() - -# Then, we define the ODE solver. That is, the scheme that will be used for the time integration. In this tutorial we use the `ThetaMethod` with $\theta = 0.5$, resulting in a 2nd order scheme. The `ThetaMethod` function receives the linear solver, the time step size $\Delta t$ (constant) and the value of $\theta$. -Δt = 0.05 -θ = 0.5 -ode_solver = ThetaMethod(linear_solver,Δt,θ) - -# Finally, we define the solution using the `solve` function, giving the ODE solver, the FE operator, an initial solution, an initial time and a final time. To construct the initial condition we interpolate the initial value (in that case a constant value of 0.0) into the FE space $U(t)$ at $t=0.0$. -u₀ = interpolate_everywhere(0.0,U(0.0)) -t₀ = 0.0 -T = 10.0 -uₕₜ = solve(ode_solver,op,u₀,t₀,T) - -# ## Postprocessing - -# We should highlight that `uₕₜ` is just an _iterable_ function and the results at each time steps are only computed when iterating over it, i.e., lazily. We can post-process the results and generate the corresponding `vtk` files using the `createpvd` and `createvtk` functions. The former will create a `.pvd` file with the collection of `.vtu` files saved at each time step by `createvtk`. The computation of the problem solutions will be triggered in the following loop: -createpvd("poisson_transient_solution") do pvd - for (uₕ,t) in uₕₜ - pvd[t] = createvtk(Ω,"poisson_transient_solution_$t"*".vtu",cellfields=["u"=>uₕ]) - end -end - -# ![](../assets/poisson_transient/poisson_transient.gif) diff --git a/src/transient_linear.jl b/src/transient_linear.jl new file mode 100644 index 00000000..21c512b2 --- /dev/null +++ b/src/transient_linear.jl @@ -0,0 +1,122 @@ +# In this tutorial, we will learn +# * How to solve a simple time-dependent PDE in Gridap +# * How to build a transient FE space +# * How to define a transient weak form +# * How to set up a time-marching scheme for a linear ODE +# * How to visualise transient results + +# We split the presentation of the ODE module of Gridap in two parts: +# * In this tutorial we focus on the differences between the steady and time-dependent cases, on a simple linear time-dependent PDE. +# * [Tutorial 18](@ref transient_nonlinear.jl) will introduce more advanced features of the ODE module of Gridap, applied to a nonlinear time-dependent PDE. + +# The [documentation of the ODE module of Gridap](https://gridap.github.io/Gridap.jl/dev/ODEs/) contains a detailed description of the framework for transient problems implemented in Gridap, including a [classification of transient problem](https://gridap.github.io/Gridap.jl/dev/ODEs/#Classification-of-ODEs), a [classification of numerical schemes](https://gridap.github.io/Gridap.jl/dev/ODEs/#Classification-of-numerical-schemes), an overview of the [high-level API](https://gridap.github.io/Gridap.jl/dev/ODEs/#High-level-API-in-Gridap) which this tutorial illustrates and of the [internals of the ODE module](https://gridap.github.io/Gridap.jl/dev/ODEs/#Low-level-implementation), as well as some [notes on and analysis of the numerical schemes implemented in Gridap](https://gridap.github.io/Gridap.jl/dev/ODEs/#Numerical-schemes-formulation-and-implementation). + +# ## Problem statement + +# We solve the heat equation in the two-dimensional domain $\Omega = [-1, +1]^{2}$. Let $k$ denote the thermal conductivity of the material, $c$ its specific heat capacity, $\rho$ its density and $q$ a rate of external heat generation. Let $g$ denote the temperature on the boundary of the domain. Let $t_{0}$ be the initial time, and $u_{0}$ be the initial temperature. The strong form of the heat equation reads: find $u(t): \Omega \to \mathbb{R}$ such that +# ```math +# \left\lbrace +# \begin{aligned} +# \rho(t, x) c(t, x) \partial_{t} u(t, x) - \nabla \cdot (k(t, x) \nabla u(t, x)) &= q(t, x) & \text{ in } \Omega, \\ +# u(t, x) &= g(t, x) & \text{ on } \partial \Omega, \\ +# u(t_{0}, x) &= u_{0}(x) & \text{ in } \Omega \\ +# \end{aligned} +# \right. +# ``` +# We assume that the data ($k$, $c$, $\rho$ and $q$) is continuous in time. Let $\alpha = k / (\rho c)$ denote the thermal diffusivity and $f = q / (\rho c)$ be the rate of external temperature generation. The weak form of the problem reads: find $u(t) \in U_{g}(t)$ such that $b(t, u, v) = \ell(t, v)$ for all $t \geq t_{0}$ and $v \in V_{0}$, where the time-dependent bilinear and linear forms $b(t, \cdot, \cdot)$ and $\ell(t, \cdot)$ are defined as +# ```math +# \begin{aligned} +# b(t, u, v) &= m(t, u, v) + a(t, u, v), \\ +# m(t, u, v) &= \int_{\Omega} v \partial_{t} u(t) \ {\rm d} \Omega, \\ +# a(t, u, v) &= \int_{\Omega} \alpha(t) \nabla v \cdot \nabla u(t) \ {\rm d} \Omega, \\ +# \ell(t, v) &= \int_{\Omega} v f(t) \ {\rm d} \Omega, +# \end{aligned} +# ``` +# and the the functional spaces are $U_{g}(t) = H^{1}_{g(t)}(\Omega)$ and $V_{0} = H^{1}_{0}(\Omega)$. In particular, the trial space $U_{g}$ is a transient functional space, in the sense that its Dirichlet trace $g$ depends on time. However, the test space $V_{0}$ is time-independent (it has a constant, zero Dirichlet trace). For all $t \geq t_{0}$, we assume that $\alpha(t) \in L^{\infty}(\Omega)$ is uniformly positive in $\Omega$, $f(t) \in H^{-1}(\Omega)$, $g(t) \in H^{1/2}(\Omega)$, and finally $u_{0} \in L^{2}(\Omega)$. + +# ## Discrete model + +# We start with the discretization of the computational domain. In our case, we consider a $20 \times 20$ Cartesian mesh of the square $[-1, +1]^{2}$. + +using Gridap +domain = (-1, +1, -1, +1) +partition = (20, 20) +model = CartesianDiscreteModel(domain, partition) + +# ## FE spaces + +# In this tutorial we use a linear Lagrangian FE space. + +order = 1 +reffe = ReferenceFE(lagrangian, Float64, order) + +# The test space is defined as for steady problems + +V0 = TestFESpace(model, reffe, dirichlet_tags="boundary") + +# The trial space is now a `TransientTrialFESpace`, which is constructed from a `TestFESpace` and a function (or vector of functions if several Dirichlet tags are provided) for the Dirichlet boundary conditions. The Dirichlet trace has to be prescribed as a function of time and then space as follows + +g(t) = x -> exp(-2 * t) * sinpi(t * x[1]) * (x[2]^2 - 1) +Ug = TransientTrialFESpace(V0, g) + +# ## Triangulation and quadrature + +# As usual, we equip the model with an integration mesh and a measure + +degree = 2 +Ω = Triangulation(model) +dΩ = Measure(Ω, degree) + +# ## Weak form +# We define the thermal diffusivity $\alpha$ and the rate of external temperature generation $f$. + +α(t) = x -> 1 + sin(t) * (x[1]^2 + x[2]^2) / 4 +f(t) = x -> sin(t) * sinpi(x[1]) * sinpi(x[2]) + +# We are going to construct a transient linear FEOperator by providing the bilinear forms associated to $\partial_{t} u$ and $u$, as well as the forcing term. Note that they now receive time as an additional argument, and the time derivative operator is `∂t`. + +m(t, dtu, v) = ∫(v * dtu)dΩ +a(t, u, v) = ∫(α(t) * ∇(v) ⋅ ∇(u))dΩ +l(t, v) = ∫(v * f(t))dΩ +op = TransientLinearFEOperator((m, a), l, Ug, V0) + +# In our case, the mass term ($m(t, \cdot, \cdot)$) is constant in time. We can take advantage of that to save some computational effort, and indicate it to Gridap as follows +op_opt = TransientLinearFEOperator((m, a), l, Ug, V0, constant_forms=(true, false)) + +# If the stiffness term ($a(t, \cdot, \cdot)$) had been constant in time, we could have set `constant_forms=(true, true)`. + +# ## Transient solver + +# Once we have defined the FE operator, we proceed with the definition of the ODE solver, i.e. the scheme that will be used for the integration in time. In this tutorial, we use the `ThetaMethod` with $\theta = 1/2$, resulting in a second-order scheme. The `ThetaMethod` function receives a solver for systems of equations, the time step size $\Delta t$ (constant) and the value of $\theta \in [0, 1]$. Since the ODE is linear the systems of equation that will arise in the time-marching scheme will be linear so we can provide `ThetaMethod` with a linear solver. + +ls = LUSolver() +Δt = 0.05 +θ = 0.5 +solver = ThetaMethod(ls, Δt, θ) + +# Gridap also implements explicit and diagonally-implicit Runge-Kutta schemes. One can access the full list of available Butcher tableaus through the exported constant `available_tableaus`. There are also constructors for explicit 2- and 3-stage schemes: `EXRK22(α)` and `EXRK33(α, β)`, `EXRK33_1(α)`, `EXRK33_2(α)` respectively, and diagonally-implicit 1- and 2-stage schemes: `SDIRK11(α)`, `SDIRK12()`, `SDIRK22(α, β, γ)`, `SDIRK23(λ)`. See the documentation of [Runge-Kutta schemes in Gridap](https://gridap.github.io/Gridap.jl/dev/ODEs/#Runge-Kutta) for a description of the corresponding tableaus. For example, one could have chosen a two-stage singly-diagonally-implicit scheme (of order 2) as follows. +tableau = :SDIRK_2_2 +solver_rk = RungeKutta(ls, ls, Δt, tableau) + +# Let $t_{F} > t_{0}$ be a final time, until when we want to evolve the problem. We define the solution using the `solve` function, giving the ODE solver, the transient FE operator, the initial and final times, and the initial solution. To construct the initial condition we interpolate the initial function $u_{0}$ onto the FE space $U_{g}$ at the initial time. In our case, $u_{0}$ is simply $g(t_{0})$. + +t0, tF = 0.0, 10.0 +uh0 = interpolate_everywhere(g(t0), Ug(t0)) +uh = solve(solver, op, t0, tF, uh0) + +# ## Postprocessing + +# We highlight that `uh` is an iterable function and the result at each time step is only computed lazily when iterating over it. We can post-process the results and generate the corresponding `vtk` files using the `createpvd` and `createvtk` functions. The former will create a `.pvd` file with the collection of `.vtu` files saved at each time step by `createvtk`. The computation of the problem solutions will be triggered in the following loop: + +if !isdir("tmp") + mkdir("tmp") +end + +createpvd("results") do pvd + pvd[0] = createvtk(Ω, "tmp/results_0" * ".vtu", cellfields=["u" => uh0]) + for (tn, uhn) in uh + pvd[tn] = createvtk(Ω, "tmp/results_$tn" * ".vtu", cellfields=["u" => uhn]) + end +end + +# ![](../assets/transient_linear/result.gif) diff --git a/src/transient_nonlinear.jl b/src/transient_nonlinear.jl new file mode 100644 index 00000000..f7baeb21 --- /dev/null +++ b/src/transient_nonlinear.jl @@ -0,0 +1,140 @@ +# In this tutorial, we will learn +# * How to write a nonlinear transient weak form in Gridap +# * How to setup a time-marching scheme for a nonlinear ODE + +# We assume that the reader is familiar with Gridap's API for linear transient PDEs, introduced in [Tutorial 17](@ref transient_linear.jl). We focus here on more advanced features of the ODE module of Gridap, applied to a nonlinear time-dependent PDE. + +# ## Problem statement + +# We consider the same problem as in [Tutorial 17](@ref transient_linear.jl), and use the same notations: find $u(t): \Omega \to \mathbb{R}$ such that +# ```math +# \left\lbrace +# \begin{aligned} +# \rho(t, x) c(t, x) \partial_{t} u(t, x) - \nabla \cdot (k(t, x) \nabla u(t, x)) &= q(t, x) & \text{ in } \Omega, \\ +# u(t, x) &= g(t, x) & \text{ on } \partial \Omega, \\ +# u(t_{0}, x) &= u_{0}(x) & \text{ in } \Omega \\ +# \end{aligned} +# \right. +# ``` +# In this tutorial we consider a nonlinear (quadratic) conductivity coefficient $\alpha(t, x, u) = \alpha_{0}(t, x) + \alpha_{1}(t, x) u(t, x) + \alpha_{2}(t, x) u(t, x)^{2}$. Here again, we assume that the $\alpha_{i}$ are continuous in time. The weak form of the problem reads: find $u(t) \in U_{g}(t)$ such that $b(t, u, v) = \ell(t, v)$ for all $t \geq t_{0}$ and $v \in V_{0}$, where the time-dependent bilinear and linear forms $b(t, \cdot, \cdot)$ and $\ell(t, \cdot)$ are defined as +# ```math +# \begin{aligned} +# b(t, u, v) &= m(t, u, v) + a(t, u, v), \\ +# m(t, u, v) &= \int_{\Omega} v \partial_{t} u(t) \ {\rm d} \Omega, \\ +# a(t, u, v) &= \int_{\Omega} \nabla v \cdot [(\alpha_{0}(t) + \alpha_{1}(t) u(t) + \alpha_{2}(t) u(t)^{2}) \nabla u(t)] \ {\rm d} \Omega, \\ +# \ell(t, v) &= \int_{\Omega} v f(t) \ {\rm d} \Omega, +# \end{aligned} +# ``` +# and the the functional spaces are $U_{g}(t) = \{u \in H^{1}_{g(t)}(\Omega), u \nabla u \in \boldsymbol{L}^{2}(\Omega), u^{2} \nabla u \in \boldsymbol{L}^{2}(\Omega)\}$ and $V_{0} = H^{1}_{0}(\Omega)$. In addition to the regularity conditions of Tutorial 17 on $f$, $g$ and $u_{0}$, we assume that for all $t \geq t_{0}$, it holds $\alpha_{i}(t) \in L^{\infty}(\Omega)$ and $(x, X) \mapsto \alpha_{0}(t, x) + \alpha_{1}(t, x) X + \alpha_{2}(t, x) X^{2}$ is uniformly positive in $\Omega \times \mathbb{R}$, i.e. $\alpha_{2}(t)$ and $4 \alpha_{0}(t) \alpha_{2}(t) - \alpha_{1}^{2}(t)$ are uniformly positive. + +# ## Discrete model, FE spaces, triangulation and quadrature + +# We consider the same mesh, FE spaces, triangulation and quadrature as in Tutorial 17: + +using Gridap +domain = (-1, +1, -1, +1) +partition = (20, 20) +model = CartesianDiscreteModel(domain, partition) + +order = 1 +reffe = ReferenceFE(lagrangian, Float64, order) + +V0 = TestFESpace(model, reffe, dirichlet_tags="boundary") + +g(t) = x -> exp(-2 * t) * sinpi(t * x[1]) * (x[2]^2 - 1) +Ug = TransientTrialFESpace(V0, g) + +degree = 2 +Ω = Triangulation(model) +dΩ = Measure(Ω, degree) + +# ## Nonlinear weak form +# We define the diffusion coefficients $\alpha$ and $\beta$, the total nonlinear diffusion coefficient $\kappa$ as well as the forcing term $f$. + +α₀(t) = x -> 1 + sin(t) * (x[1]^2 + x[2]^2) / 4 +α₁(t) = x -> cos(t) * x[1]^2 / 2 +α₂(t) = x -> 1 + t * (x[1]^2 + x[2]^2) +α(t, u) = α₀(t) + α₁(t) * u + α₂(t) * u * u +f(t) = x -> sin(t) * sinpi(x[1]) * sinpi(x[2]) + +# We now write the nonlinear weak form. Similar to steady nonlinear problems, we provide the residual and its Jacobian, here with respect to $u$ and $\partial_{t} u$. The mass, stiffness and forcing terms are written as follows. + +m(t, u, v) = ∫(v * ∂t(u))dΩ +a(t, u, v) = ∫(∇(v) ⋅ (α(t, u) * ∇(u)))dΩ +l(t, v) = ∫(v * f(t))dΩ + +# The Jacobians of the mass and the stiffness are + +jac_m(t, u, dtu, v) = ∫(v * dtu)dΩ +jac_α(t, u, du) = α₁(t) * du + α₂(t) * (2 * u * du) +jac_a(t, u, du, v) = ∫(∇(v) ⋅ (α(t, u) * ∇(du)))dΩ + ∫(∇(v) ⋅ (jac_α(t, u, du) * ∇(u)))dΩ + +# We can now write the residual and its Jacobians with respect to $u$ and $\partial_{t} u$ as follows + +res(t, u, v) = m(t, u, v) + a(t, u, v) - l(t, v) +jac(t, u, du, v) = jac_a(t, u, du, v) +jac_t(t, u, dtu, v) = jac_m(t, u, dtu, v) + +# The most general way of constructing a transient FE operator is by using the `TransientFEOperator` constructor, which receives a residual, a Jacobian with respect to the unknown and a Jacobian with respect to the time derivative. + +op = TransientFEOperator(res, (jac, jac_t), Ug, V0) + +# In this example, the mass term is linear so this ODE belongs to the class of quasilinear ODEs. We can indicate this additional structure to Gridap as follows + +mass_ql(t, u, dtu, v) = ∫(dtu * v)dΩ +res_ql(t, u, v) = a(t, u, v) - l(t, v) +jac_ql(t, u, du, v) = jac_a(t, u, du, v) +jac_t_ql(t, u, dtu, v) = jac_m(t, u, dtu, v) +op_ql = TransientQuasilinearFEOperator(mass_ql, res_ql, (jac_ql, jac_t_ql), Ug, V0) + +# In fact, this ODE further classifies as semilinear because its mass term does not involve $u$. In our case, the mass term is also constant in time, so the optimal operator is as follows. Note that the signature of the mass term does not involve `u` anymore, as this is the condition for an ODE to be semilinear. + +mass_sl(t, dtu, v) = ∫(dtu * v)dΩ +res_sl(t, u, v) = a(t, u, v) - l(t, v) +jac_sl(t, u, du, v) = jac_a(t, u, du, v) +jac_t_sl(t, u, dtu, v) = mass_sl(t, dtu, v) +op_sl = TransientSemilinearFEOperator( + mass_sl, res_sl, (jac_sl, jac_t_sl), + Ug, V0, constant_mass=true +) + +# In all cases above, it is also possible to take advantage of automatic differentiation techniques to compute both Jacobians and build the transient FE operator from the residual and the FE spaces only. + +# ## Transient solver + +# We proceed to the definition of the ODE solver. If the ODE is described via a general nonlinear FE operator, we will need to provide these schemes with a nonlinear solver for systems of equations. If the operator is quasilinear and the scheme is explicit, one only needs a linear solver. Here we draw from `NLSolvers.jl` and rely on a Newton-Raphson solver based on Gridap's `LUSolver`. + +# For example, for the `ThetaMethod`, one would write +lin_solver = LUSolver() +nl_solver = NLSolver(lin_solver, method=:newton, iterations=10, show_trace=false) + +Δt = 0.05 +θ = 0.5 +solver = ThetaMethod(nl_solver, Δt, θ) + +# For a two-stage singly-diagonally-implicit scheme (of order 2), it would be +tableau = :SDIRK_2_2 +solver_rk = RungeKutta(nl_solver, lin_solver, Δt, tableau) + +# We define the initial condition and the solution using the `solve` function as in Tutorial 17: + +t0, tF = 0.0, 10.0 +uh0 = interpolate_everywhere(g(t0), Ug(t0)) +uh = solve(solver, op_sl, t0, tF, uh0) + +# ## Postprocessing + +# Here again, we export the solution at each time step as follows + +if !isdir("tmp_nl") + mkdir("tmp_nl") +end + +createpvd("results_nl") do pvd + pvd[0] = createvtk(Ω, "tmp_nl/results_0" * ".vtu", cellfields=["u" => uh0]) + for (tn, uhn) in uh + pvd[tn] = createvtk(Ω, "tmp_nl/results_$tn" * ".vtu", cellfields=["u" => uhn]) + end +end + +# ![](../assets/transient_nonlinear/result.gif)