# Method of Manufactured Solutions for non-differentiable coefficients

The Method of manufactured is commonly used for verification of computational science codes. I want to use the method for verification of Navier-Cauchy (elasticity) equations with periodic and Bloch-periodic boundary conditions. In this case, the coefficients of the PDE (material propertise) are periodic, and the form of the boundary value problem is

$$\nabla\cdot \sigma + \mathbf{f} = -\omega^2 \rho \mathbf{u}\quad \forall \mathbf{x} \in \Omega\\ \mathbf{u}(\mathbf{x} + \mathbf{a}) = \mathbf{u}(\mathbf{x})e^{\mathbf{k}\cdot\mathbf{a}}\quad \forall \mathbf{x} \in \Gamma_u\\ \mathbf{t}(\mathbf{x} + \mathbf{a}) = -\mathbf{t}(\mathbf{x})e^{\mathbf{k}\cdot\mathbf{a}}\quad \forall \mathbf{x} \in \Gamma_t$$ with \begin{align} &\sigma = \underline{C} \epsilon\\ &\epsilon = \frac{1}{2}\left[\nabla\mathbf{u} + \nabla\mathbf{u}^T\right] \enspace , \end{align} and $\underline{C}$ presents the desired periodicity, i.e., $\underline{C}(\mathbf{x} + \mathbf{a}) = \underline{C}(\mathbf{x})$, being $\mathbf{a}$ the vector periodicity of the material.

When $\omega=0$ the Bloch-periodic conditions turn into common periodic conditions and $\mathbf{u}(\mathbf{x} + \mathbf{a}) = \mathbf{u}(\mathbf{x})$ (similar for the tractions $\mathbf{t}(\mathbf{x})$). In that case one can propose a solution that presents the same periodicity of the material, e.g

$$\mathbf{u} = \left[e^{\sin \left(2\,x\right)\,\sin y}-\cos \left(3\,x\right)\,\sin ^2 \left(2\,y\right)\right] (u_1, u_2, 0)$$ and the periodicity for the material being $$\underline{C} = \left[\sin ^2x\,\sin \cos y^3+\sin \cos x^2\,\sin ^4y\right] \underline{C}_0 \enspace ,$$

then, plug these terms in the equations and obtain the expression for $\mathbf{f}$ (this can be done symbolically using Maxima or SymPy and then generate Fortran, C++ or Python code).

In the case of Bloch-periodicity it is a little bit more convoluted, since the functions that solve the equation are not periodic. They are the product of two periodic functions, and just periodic in the particular case of commensurate wavelengths. Thus, in this case the solution would depend on the wavevector $\mathbf{k}$.

Question: It is customary to assume that the material properties are composites where the properties change abruptly, e.g., a bilayer material, fiber reinforced composites. In these cases the constitutive tensor is not a differentiable function. If the code is verified with smooth functions what can be said about the case of non-differentiable material properties? Is the test still valid?

• I think you can do you MMS through the weak form if you don't try to demand a classical solution. You still end up with some terms to implement for the RHS. – Bill Barth Jul 6 '15 at 21:12
• @BillBarth That sounds interesting, but not easy. I just tried to compute the RHS for the same example proposed up and constant coefficients and the integral cannot be solved analytically in Maxima or Maple. This is the simplest case, so I think that the approach cannot be applied (easily) for more general cases. – nicoguaro Jul 6 '15 at 22:41
• For your second question, the test may be valid for smooth solutions (though you will likely need to implement the RHS with Bill's technique, perhaps numerically using quadrature). For example, with Poisson with discontinuous coeffs, FEM should still achieve optimal rates of convergence if the discontinuities are fitted to the mesh. Is this what you mean by valid? – Jesse Chan Jul 6 '15 at 23:16
• @JesseChan I was referring to compute the verification for smooth cases and expect that the code will work correctly for the non-smooth case (since the functions are piecewise continuous and differentiable). Regarding the numerical integration, is there a way to achieve the quadrature for different meshes in an accurate way? – nicoguaro Jul 6 '15 at 23:23
• Since integration is typically done locally, I would expect that to work (i.e. there is no assumption of smoothness in computations of integrals). For quadrature, I only know about finite elements and convergence tests. Cubature points are mapped from a reference element to the physical domain, and it usually suffices to take a cubature rule of one degree higher, since cubature errors decrease at a faster rate relative to expected discretization errors. – Jesse Chan Jul 7 '15 at 0:11

Let me see if I can turn my comment into a fully-fledged answer. Suppose you are solving the Laplacian with a variable $k$ which has a step function across the domain. Suppose you have homogeneous boundary conditions so that we don't have to deal with that in my answer. Now your weak form will look something like:

$$\int_\Omega k\nabla u \cdot \nabla v \; dx = \int_\Omega f v \;dx \;\;\forall v$$

Suppose you want to have $u=\tanh(x)\sin(x)$, in the strong form, you can't differentiate $k$, because it's discontinuous, but you can compute

$$\int_\Omega k\nabla (\tanh(x)\sin(x)) \cdot \nabla v \; dx = \int_\Omega f v \;dx \;\;\forall v$$

Now, given that, instead of trying to compute $f$ from this, you can just implement the above as your right-hand side. I.e., solve

$$\int_\Omega k\nabla u_h \cdot \nabla v_h \; dx = \int_\Omega k\nabla (\tanh(x)\sin(x)) \cdot \nabla v_h \; dx \;\;\forall v_h$$

using your favorite $k$. As long as your integration makes sense, and as long as you can take that single derivative of your chosen solution, you should be good to go.

• Thanks! Then, I think that the way to go is to perform both the smooth case with complicated coefficients and another one with the "composite case" where the properties are not differentiable. And, as far as I understand from your comments and answer, in any case it is necessary to compute some of the integrals using quadrature. – nicoguaro Jul 8 '15 at 18:25
• @nicoguaro, I'm not sure I understand your response. The thing with MMS is that you actually compute something so of course you're going to do quadrature at some point. In any complex PDE, you're not going to be able to take an arbitrary assumed solution and exactly integrate the necessary RHS. Quadrature is a part of the process. – Bill Barth Jul 8 '15 at 18:32
• On one hand, I was asking about the quadrature, and you already reply this in the comment (I read some papers and reports about the MMS, and they explained things for FDM... so they never mentioned that part). On the other hand, I was mentioning that is worth to consider one smooth (and convoluted) function for the coefficients and also a non-differentiable (but simpler) one to get the "best" test. – nicoguaro Jul 8 '15 at 18:36