diff --git a/src/TopOptEMFocus.jl b/src/TopOptEMFocus.jl index b526d25a..9933e0de 100644 --- a/src/TopOptEMFocus.jl +++ b/src/TopOptEMFocus.jl @@ -70,7 +70,7 @@ # ``` # 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 lengthscale $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$: +# 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} # -r_f^2\nabla^2p_f+p_f&=p\, ,\\ @@ -182,14 +182,14 @@ 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. # @@ -318,7 +318,7 @@ end # ## Optimization with adjoint method # -# Now that we have our objective to optimize, the next step is to find out the derivative to the desigin 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 diff --git a/src/darcy.jl b/src/darcy.jl index 420e79af..420eb4df 100644 --- a/src/darcy.jl +++ b/src/darcy.jl @@ -85,13 +85,13 @@ X = MultiFieldFESpace([U, P]) # ## Numerical integration # -# In this example we need to integrate in the interior of $\Omega$ and on the Neumann boundary $\Gamma_{\rm N}$. For the volume integrals, we extract the triangulation from the geometrical model and define the corresponding Lebesge measures, which will allow to write down the integrals of the weak form. +# In this example we need to integrate in the interior of $\Omega$ and on the Neumann boundary $\Gamma_{\rm N}$. For the volume integrals, we extract the triangulation from the geometrical model and define the corresponding Lebesgue measures, which will allow to write down the integrals of the weak form. trian = Triangulation(model) degree = 2 dΩ = Measure(trian,degree) -# In order to integrate the Neumann boundary condition, we only need to build an integration mesh for the right side of the domain (which is the only part of $\Gamma_{\rm N}$, where the Neumann function $h$ is different from zero). Within the model, the right side of $\Omega$ is identified with the boundary tag 8. Using this identifier, we extract the corresponding surface triangulation and create the required Lebesge measure. +# In order to integrate the Neumann boundary condition, we only need to build an integration mesh for the right side of the domain (which is the only part of $\Gamma_{\rm N}$, where the Neumann function $h$ is different from zero). Within the model, the right side of $\Omega$ is identified with the boundary tag 8. Using this identifier, we extract the corresponding surface triangulation and create the required Lebesgue measure. neumanntags = [8,] btrian = BoundaryTriangulation(model,tags=neumanntags) diff --git a/src/dg_discretization.jl b/src/dg_discretization.jl index d309d181..7b444a74 100644 --- a/src/dg_discretization.jl +++ b/src/dg_discretization.jl @@ -245,7 +245,7 @@ writevtk(Λ,"strian") # ![](../assets/dg_discretization/skeleton_trian.png) # # Once we have constructed the triangulations needed in this example, we define -# the corresponding quadrature rules by passing the triangualtions +# the corresponding quadrature rules by passing the triangulations # together with the desired degree to the `Measure` function. degree = 2*order diff --git a/src/emscatter.jl b/src/emscatter.jl index 14034b14..570ee3ba 100644 --- a/src/emscatter.jl +++ b/src/emscatter.jl @@ -35,7 +35,7 @@ # 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 in $x/y$ direction, the associated boundary condition is then usually the zero Dirichlet boundary condition. Here, the boundary conditions are zero Dirichlet boundary on the top and bottom side $\Gamma_D$ but periodic boundary condition on the left ($\Gamma_L$) and right side ($\Gamma_R$). -# The reason that we use a periodic boundary condition for the left and right side instead of zero Dirichlet boundary condition is that we want to simulate a plane wave exicitation, which then requires a periodic boundary condition. +# The reason that we use a periodic boundary condition for the left and right side instead of zero Dirichlet boundary condition is that we want to simulate a plane wave excitation, which then requires a periodic boundary condition. # # Consider $\mu(x)=1$ (which is mostly the case in electromagnetic problems) 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 formulate the problem as # @@ -46,7 +46,7 @@ # 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 instead of a diagonal $2 \times 2$ matrix, in which case $\Lambda\nabla$ becomes the elementwise product. +# For convenience, in the weak form and Julia implementation below we represent $\Lambda$ as a vector instead of a diagonal $2 \times 2$ matrix, in which case $\Lambda\nabla$ becomes the element-wise product. # # ## Numerical scheme # @@ -89,7 +89,7 @@ model = GmshDiscreteModel("../models/geometry.msh") # ## FE spaces # -# We use the first-order lagrangian as the finite element function space basis. The Dirichlet edges are labeled with `DirichletEdges` in the mesh file. Since our problem involves complex numbers (because of PML), we need to assign the `vector_type` to be `Vector{ComplexF64}`. +# We use the first-order Lagrangian as the finite element function space basis. The Dirichlet edges are labeled with `DirichletEdges` in the mesh file. Since our problem involves complex numbers (because of PML), we need to assign the `vector_type` to be `Vector{ComplexF64}`. # # ### Test and trial finite element function space @@ -120,7 +120,7 @@ dΓ = Measure(Γ,degree) # ### PML parameters Rpml = 1e-12 # Tolerance for PML reflection σ = -3/4*log(Rpml)/d_pml # σ_0 -LH = (L,H) # Size of the PML inner boundary (a rectangular centere at (0,0)) +LH = (L,H) # Size of the PML inner boundary (a rectangular center at (0,0)) # ### PML coordinate stretching functions function s_PML(x,σ,k,LH,d_pml) @@ -189,7 +189,7 @@ uh = solve(op) # ## Analytical solution # ### Theoretical analysis -# In this section, we construct the semi-analytical solution to this scattering problem, for comparison to the numerical solution. This is possible because of the symmetry of the cylinder, which allows us to expand the solutions of the Helmoltz equation in Bessel functions and match boundary conditions at the cylinder interface. (In 3d, the analogous process with spherical harmonics is known as "Mie scattering".) For more information on this technique, see Ref [4]. +# In this section, we construct the semi-analytical solution to this scattering problem, for comparison to the numerical solution. This is possible because of the symmetry of the cylinder, which allows us to expand the solutions of the Helmholtz equation in Bessel functions and match boundary conditions at the cylinder interface. (In 3d, the analogous process with spherical harmonics is known as "Mie scattering".) For more information on this technique, see Ref [4]. # In 2D cylinder coordinates, we can expand the plane wave in terms of Bessel functions (this is the Jacobi–Anger identity [5]): # # ```math diff --git a/src/fsi_tutorial.jl b/src/fsi_tutorial.jl index 8b4b722b..cc996ae4 100644 --- a/src/fsi_tutorial.jl +++ b/src/fsi_tutorial.jl @@ -189,7 +189,7 @@ n_Γout = get_normal_vector(Γ_out) Γ_fs = InterfaceTriangulation(Ω_f,Ω_s) n_Γfs = get_normal_vector(Γ_fs) -# In what follows we will assume a second-order veloticty interpolation, i.e. $k=2$ +# In what follows we will assume a second-order velocity interpolation, i.e. $k=2$ k = 2 # Now we define the reference FE for the velocity and pressure fields. The velocity field reference FE, both for fluid and solid domains, will be defined by a 2-dimensional `VectorValue` type `:Lagrangian` reference FE element of order `k`. The reference FE for the pressure will be given by a scalar value type `:Lagrangian` reference FE element of order `k-1`. diff --git a/src/poisson_transient.jl b/src/poisson_transient.jl index e53c58f9..7a473d3c 100644 --- a/src/poisson_transient.jl +++ b/src/poisson_transient.jl @@ -63,7 +63,7 @@ U = TransientTrialFESpace(V,g) # 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. +# 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Ω diff --git a/src/validation.jl b/src/validation.jl index 31041e02..5c52be93 100644 --- a/src/validation.jl +++ b/src/validation.jl @@ -107,7 +107,7 @@ uh = solve(op) # # ## Measuring the discretization error # -# Our goal is to check that the discratization error associated with the computed approximation `uh` is close to machine precision. To this end, the first step is to compute the discretization error, which is done as you would expect: +# Our goal is to check that the discretization error associated with the computed approximation `uh` is close to machine precision. To this end, the first step is to compute the discretization error, which is done as you would expect: e = u - uh diff --git a/src/validation_DrWatson.jl b/src/validation_DrWatson.jl index 0b3916ac..92a4b40e 100644 --- a/src/validation_DrWatson.jl +++ b/src/validation_DrWatson.jl @@ -1,7 +1,7 @@ # In this tutorial, we will learn # - How to use `DrWatson.jl` to accelerate and reproduce our Gridap simulation workflows -# [DrWatson.jl](https://github.com/JuliaDynamics/DrWatson.jl) is a Julia package that helps managing a typical scientific workflow thorought all its phases, see a summary [here](https://juliadynamics.github.io/DrWatson.jl/stable/workflow/). +# [DrWatson.jl](https://github.com/JuliaDynamics/DrWatson.jl) is a Julia package that helps managing a typical scientific workflow through all its phases, see a summary [here](https://juliadynamics.github.io/DrWatson.jl/stable/workflow/). # All its functionalities can be accessed with (non-invasive) simple function calls. @@ -40,7 +40,7 @@ import Gridap: ∇ # To this end, we want to solve our computational model for many combinations of mesh size and order of FE approximation (*parameters*) and extract the L2- and H1-norm errors (*output data*). -# We first group all parameters and parameter values in a single ditionary +# We first group all parameters and parameter values in a single dictionary params = Dict( "cells_per_axis" => [8,16,32,64],