ex13.php

一个用来实现偏微分方程中网格的计算库

  
  
  
  <B><FONT COLOR="#228B22">void</FONT></B> assemble_stokes (EquationSystems&amp; es,
  		      <B><FONT COLOR="#228B22">const</FONT></B> std::string&amp; system_name)
  {
    assert (system_name == <B><FONT COLOR="#BC8F8F">&quot;Navier-Stokes&quot;</FONT></B>);
    
    <B><FONT COLOR="#228B22">const</FONT></B> Mesh&amp; mesh = es.get_mesh();
    
    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dim = mesh.mesh_dimension();
    
    TransientLinearImplicitSystem &amp; navier_stokes_system =
      es.get_system&lt;TransientLinearImplicitSystem&gt; (<B><FONT COLOR="#BC8F8F">&quot;Navier-Stokes&quot;</FONT></B>);
  
    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> u_var = navier_stokes_system.variable_number (<B><FONT COLOR="#BC8F8F">&quot;u&quot;</FONT></B>);
    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> v_var = navier_stokes_system.variable_number (<B><FONT COLOR="#BC8F8F">&quot;v&quot;</FONT></B>);
    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> p_var = navier_stokes_system.variable_number (<B><FONT COLOR="#BC8F8F">&quot;p&quot;</FONT></B>);
    
    FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
    
    FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
  
    AutoPtr&lt;FEBase&gt; fe_vel  (FEBase::build(dim, fe_vel_type));
      
    AutoPtr&lt;FEBase&gt; fe_pres (FEBase::build(dim, fe_pres_type));
    
    QGauss qrule (dim, fe_vel_type.default_quadrature_order());
  
    fe_vel-&gt;attach_quadrature_rule (&amp;qrule);
    fe_pres-&gt;attach_quadrature_rule (&amp;qrule);
    
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;Real&gt;&amp; JxW = fe_vel-&gt;get_JxW();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;Real&gt; &gt;&amp; phi = fe_vel-&gt;get_phi();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;RealGradient&gt; &gt;&amp; dphi = fe_vel-&gt;get_dphi();
  
    <B><FONT COLOR="#228B22">const</FONT></B> std::vector&lt;std::vector&lt;Real&gt; &gt;&amp; psi = fe_pres-&gt;get_phi();
  
    
    <B><FONT COLOR="#228B22">const</FONT></B> DofMap &amp; dof_map = navier_stokes_system.get_dof_map();
  
    DenseMatrix&lt;Number&gt; Ke;
    DenseVector&lt;Number&gt; Fe;
  
    DenseSubMatrix&lt;Number&gt;
      Kuu(Ke), Kuv(Ke), Kup(Ke),
      Kvu(Ke), Kvv(Ke), Kvp(Ke),
      Kpu(Ke), Kpv(Ke), Kpp(Ke);
  
    DenseSubVector&lt;Number&gt;
      Fu(Fe),
      Fv(Fe),
      Fp(Fe);
  
    <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector&lt;<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>&gt; dof_indices;
    <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector&lt;<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>&gt; dof_indices_u;
    <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector&lt;<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>&gt; dof_indices_v;
    <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector&lt;<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>&gt; dof_indices_p;
  
    <B><FONT COLOR="#228B22">const</FONT></B> Real dt    = es.parameters.get&lt;Real&gt;(<B><FONT COLOR="#BC8F8F">&quot;dt&quot;</FONT></B>);
    <B><FONT COLOR="#228B22">const</FONT></B> Real theta = 1.;
      
    <B><FONT COLOR="#5F9EA0">MeshBase</FONT></B>::const_element_iterator       el     = mesh.active_local_elements_begin();
    <B><FONT COLOR="#228B22">const</FONT></B> MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); 
    
    <B><FONT COLOR="#A020F0">for</FONT></B> ( ; el != end_el; ++el)
      {    
        <B><FONT COLOR="#228B22">const</FONT></B> Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
        dof_map.dof_indices (elem, dof_indices_u, u_var);
        dof_map.dof_indices (elem, dof_indices_v, v_var);
        dof_map.dof_indices (elem, dof_indices_p, p_var);
  
        <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_dofs   = dof_indices.size();
        <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_u_dofs = dof_indices_u.size(); 
        <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_v_dofs = dof_indices_v.size(); 
        <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_p_dofs = dof_indices_p.size();
        
        fe_vel-&gt;reinit  (elem);
        fe_pres-&gt;reinit (elem);
  
        Ke.resize (n_dofs, n_dofs);
        Fe.resize (n_dofs);
  
        Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
        Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
        Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
        
        Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
        Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
        Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
        
        Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
        Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
        Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
  
        Fu.reposition (u_var*n_u_dofs, n_u_dofs);
        Fv.reposition (v_var*n_u_dofs, n_v_dofs);
        Fp.reposition (p_var*n_u_dofs, n_p_dofs);
  
        <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp&lt;qrule.n_points(); qp++)
  	{
  	  Number   u = 0., u_old = 0.;
  	  Number   v = 0., v_old = 0.;
  	  Number   p_old = 0.;
  	  Gradient grad_u, grad_u_old;
  	  Gradient grad_v, grad_v_old;
  	  
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> l=0; l&lt;n_u_dofs; l++)
  	    {
  	      u_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_u[l]);
  	      v_old += phi[l][qp]*navier_stokes_system.old_solution (dof_indices_v[l]);
  	      grad_u_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_u[l]));
  	      grad_v_old.add_scaled (dphi[l][qp],navier_stokes_system.old_solution (dof_indices_v[l]));
  
  	      u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]); 
  	      v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
  	      grad_u.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_u[l]));
  	      grad_v.add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_v[l]));
  	    }
  
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> l=0; l&lt;n_p_dofs; l++)
  	    {
  	      p_old += psi[l][qp]*navier_stokes_system.old_solution (dof_indices_p[l]);
  	    }
  
  	  <B><FONT COLOR="#228B22">const</FONT></B> NumberVectorValue U_old (u_old, v_old);
  	  <B><FONT COLOR="#228B22">const</FONT></B> NumberVectorValue U     (u,     v);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Number  u_x = grad_u(0);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Number  u_y = grad_u(1);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Number  v_x = grad_v(0);
  	  <B><FONT COLOR="#228B22">const</FONT></B> Number  v_y = grad_v(1);
  	  
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;n_u_dofs; i++)
  	    {
  	      Fu(i) += JxW[qp]*(u_old*phi[i][qp] -                            <I><FONT COLOR="#B22222">// mass-matrix term 
</FONT></I>  				(1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] + <I><FONT COLOR="#B22222">// convection term
</FONT></I>  				(1.-theta)*dt*p_old*dphi[i][qp](0)  -         <I><FONT COLOR="#B22222">// pressure term on rhs
</FONT></I>  				(1.-theta)*dt*(grad_u_old*dphi[i][qp]) +      <I><FONT COLOR="#B22222">// diffusion term on rhs
</FONT></I>  				theta*dt*(U*grad_u)*phi[i][qp]);              <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  
  		
  	      Fv(i) += JxW[qp]*(v_old*phi[i][qp] -                             <I><FONT COLOR="#B22222">// mass-matrix term
</FONT></I>  				(1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] +  <I><FONT COLOR="#B22222">// convection term
</FONT></I>  				(1.-theta)*dt*p_old*dphi[i][qp](1) -           <I><FONT COLOR="#B22222">// pressure term on rhs
</FONT></I>  				(1.-theta)*dt*(grad_v_old*dphi[i][qp]) +       <I><FONT COLOR="#B22222">// diffusion term on rhs
</FONT></I>  				theta*dt*(U*grad_v)*phi[i][qp]);               <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  	    
  
  
  	      <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j&lt;n_u_dofs; j++)
  		{
  		  Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                <I><FONT COLOR="#B22222">// mass matrix term
</FONT></I>  				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   <I><FONT COLOR="#B22222">// diffusion term
</FONT></I>  				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  <I><FONT COLOR="#B22222">// convection term
</FONT></I>  				       theta*dt*u_x*phi[i][qp]*phi[j][qp]);   <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  
  		  Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp];     <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  		  
  		  Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +                <I><FONT COLOR="#B22222">// mass matrix term
</FONT></I>  				       theta*dt*(dphi[i][qp]*dphi[j][qp]) +   <I><FONT COLOR="#B22222">// diffusion term
</FONT></I>  				       theta*dt*(U*dphi[j][qp])*phi[i][qp] +  <I><FONT COLOR="#B22222">// convection term
</FONT></I>  				       theta*dt*v_y*phi[i][qp]*phi[j][qp]);   <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  
  		  Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp];     <I><FONT COLOR="#B22222">// Newton term
</FONT></I>  		}
  
  	      <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j&lt;n_p_dofs; j++)
  		{
  		  Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
  		  Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
  		}
  	    }
  
  	  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i&lt;n_p_dofs; i++)
  	    {
  	      <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j&lt;n_u_dofs; j++)
  		{
  		  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
  		  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
  		}
  	    }
  	} <I><FONT COLOR="#B22222">// end of the quadrature point qp-loop
</FONT></I>  
        
        {
  	<B><FONT COLOR="#228B22">const</FONT></B> Real penalty = 1.e10;
  		  
  	<B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> s=0; s&lt;elem-&gt;n_sides(); s++)
  	  <B><FONT COLOR="#A020F0">if</FONT></B> (elem-&gt;neighbor(s) == NULL)
  	    {
  	      <B><FONT COLOR="#228B22">short</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> bc_id = mesh.boundary_info-&gt;boundary_id (elem,s);
  	      <B><FONT COLOR="#A020F0">if</FONT></B> (bc_id==BoundaryInfo::invalid_id)
  		  error();
  
  	      
  	      AutoPtr&lt;Elem&gt; side (elem-&gt;build_side(s));
  	      	      
  	      <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> ns=0; ns&lt;side-&gt;n_nodes(); ns++)
  		{
  		   
  		  <B><FONT COLOR="#228B22">const</FONT></B> Real u_value = (bc_id==2) ? 1. : 0.;
  		  
  		  <B><FONT COLOR="#228B22">const</FONT></B> Real v_value = 0.;
  		  
  		  <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n=0; n&lt;elem-&gt;n_nodes(); n++)
  		    <B><FONT COLOR="#A020F0">if</FONT></B> (elem-&gt;node(n) == side-&gt;node(ns))
  		      {
  			Kuu(n,n) += penalty;
  			Kvv(n,n) += penalty;
  		  		  
  			Fu(n) += penalty*u_value;
  			Fv(n) += penalty*v_value;
  		      }
  		} <I><FONT COLOR="#B22222">// end face node loop	  
</FONT></I>  	    } <I><FONT COLOR="#B22222">// end if (elem-&gt;neighbor(side) == NULL)
</FONT></I>  	
  	<B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">bool</FONT></B> pin_pressure = true;
  	<B><FONT COLOR="#A020F0">if</FONT></B> (pin_pressure)
  	  {
  	    <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> pressure_node = 0;
  	    <B><FONT COLOR="#228B22">const</FONT></B> Real p_value               = 0.0;
  	    <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> c=0; c&lt;elem-&gt;n_nodes(); c++)
  	      <B><FONT COLOR="#A020F0">if</FONT></B> (elem-&gt;node(c) == pressure_node)
  		{
  		  Kpp(c,c) += penalty;
  		  Fp(c)    += penalty*p_value;
  		}
  	  }
        } <I><FONT COLOR="#B22222">// end boundary condition section	  
</FONT></I>        
        navier_stokes_system.matrix-&gt;add_matrix (Ke, dof_indices);
        navier_stokes_system.rhs-&gt;add_vector    (Fe, dof_indices);
      } <I><FONT COLOR="#B22222">// end of element loop
</FONT></I>    
    <B><FONT COLOR="#A020F0">return</FONT></B>;
  }
</pre> 
<a name="output"></a> 
<br><br><br> <h1> The console output of the program: </h1> 
<pre>
***************************************************************
* Running Example  ./ex13-devel
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=1681
  n_elem()=400
   n_local_elem()=400
   n_active_elem()=400
  n_subdomains()=1
  n_processors()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System "Navier-Stokes"
    Type "TransientLinearImplicit"
    Variables="u" "v" "p" 
    Finite Element Types="LAGRANGE" "LAGRANGE" "LAGRANGE" 
    Approximation Orders="SECOND" "SECOND" "FIRST" 
    n_dofs()=3803
    n_local_dofs()=3803
    n_constrained_dofs()=0
    n_vectors()=3



*** Solving time step 0, time = 0.01 ***
Linear solver converged at step: 236, final residual: 0.000214866  Nonlinear convergence: ||u - u_old|| = 282.453
Linear solver converged at step: 90, final residual: 9.99752e-06  Nonlinear convergence: ||u - u_old|| = 0.743588
Linear solver converged at step: 250, final residual: 1.04334e-07  Nonlinear convergence: ||u - u_old|| = 0.0151695
Linear solver converged at step: 250, final residual: 4.20802e-10  Nonlinear convergence: ||u - u_old|| = 0.000178465
 Nonlinear solver converged at step 3


*** Solving time step 1, time = 0.02 ***
Linear solver converged at step: 148, final residual: 0.000196443  Nonlinear convergence: ||u - u_old|| = 30.6968
Linear solver converged at step: 74, final residual: 8.13638e-06  Nonlinear convergence: ||u - u_old|| = 0.194999
Linear solver converged at step: 250, final residual: 3.57667e-08  Nonlinear convergence: ||u - u_old|| = 0.0088697
Linear solver converged at step: 250, final residual: 8.44767e-13  Nonlinear convergence: ||u - u_old|| = 5.9269e-05
 Nonlinear solver converged at step 3


*** Solving time step 2, time = 0.03 ***
Linear solver converged at step: 77, final residual: 0.000223335  Nonlinear convergence: ||u - u_old|| = 6.00665
Linear solver converged at step: 118, final residual: 9.30581e-06  Nonlinear convergence: ||u - u_old|| = 0.268524
Linear solver converged at step: 250, final residual: 3.29197e-08  Nonlinear convergence: ||u - u_old|| = 0.00316201
Linear solver converged at step: 250, final residual: 2.45713e-09  Nonlinear convergence: ||u - u_old|| = 5.82116e-05
 Nonlinear solver converged at step 3


*** Solving time step 3, time = 0.04 ***
Linear solver converged at step: 72, final residual: 0.000210495  Nonlinear convergence: ||u - u_old|| = 2.18249
Linear solver converged at step: 159, final residual: 9.538e-06  Nonlinear convergence: ||u - u_old|| = 0.260483
Linear solver converged at step: 250, final residual: 2.18784e-07  Nonlinear convergence: ||u - u_old|| = 0.0103119
Linear solver converged at step: 250, final residual: 3.4177e-09  Nonlinear convergence: ||u - u_old|| = 0.000386804
 Nonlinear solver converged at step 3


*** Solving time step 4, time = 0.05 ***
Linear solver converged at step: 58, final residual: 0.000192962  Nonlinear convergence: ||u - u_old|| = 1.02052
Linear solver converged at step: 119, final residual: 7.33031e-06  Nonlinear convergence: ||u - u_old|| = 0.103956
Linear solver converged at step: 225, final residual: 1.16591e-08  Nonlinear convergence: ||u - u_old|| = 0.00534628
Linear solver converged at step: 250, final residual: 1.41694e-12  Nonlinear convergence: ||u - u_old|| = 3.60742e-06
 Nonlinear solver converged at step 3


*** Solving time step 5, time = 0.06 ***
Linear solver converged at step: 52, final residual: 0.000214167  Nonlinear convergence: ||u - u_old|| = 0.476324
Linear solver converged at step: 172, final residual: 1.00119e-05  Nonlinear convergence: ||u - u_old|| = 0.202403
Linear solver converged at step: 158, final residual: 2.16565e-08  Nonlinear convergence: ||u - u_old|| = 0.0173926
Linear solver converged at step: 250, final residual: 2.6449e-11  Nonlinear convergence: ||u - u_old|| = 2.11892e-05
 Nonlinear solver converged at step 3


*** Solving time step 6, time = 0.07 ***
Linear solver converged at step: 21, final residual: 0.000225026  Nonlinear convergence: ||u - u_old|| = 0.256854
Linear solver converged at step: 117, final residual: 1.12946e-05  Nonlinear convergence: ||u - u_old|| = 0.147309
Linear solver converged at step: 250, final residual: 1.54091e-06  Nonlinear convergence: ||u - u_old|| = 0.0128016
Linear solver converged at step: 250, final residual: 2.11486e-09  Nonlinear convergence: ||u - u_old|| = 0.00239393
Linear solver converged at step: 250, final residual: 8.06203e-12  Nonlinear convergence: ||u - u_old|| = 2.67758e-06
 Nonlinear solver converged at step 4


*** Solving time step 7, time = 0.08 ***
Linear solver converged at step: 19, final residual: 0.000192228  Nonlinear convergence: ||u - u_old|| = 0.152223
Linear solver converged at step: 67, final residual: 7.02725e-06  Nonlinear convergence: ||u - u_old|| = 0.0850415
Linear solver converged at step: 194, final residual: 1.08592e-08  Nonlinear convergence: ||u - u_old|| = 0.00357155
Linear solver converged at step: 250, final residual: 6.3798e-12  Nonlinear convergence: ||u - u_old|| = 1.67316e-05
 Nonlinear solver converged at step 3


*** Solving time step 8, time = 0.09 ***
Linear solver converged at step: 15, final residual: 0.000204876  Nonlinear convergence: ||u - u_old|| = 0.086678
Linear solver converged at step: 97, final residual: 8.76745e-06  Nonlinear convergence: ||u - u_old|| = 0.0454467
Linear solver converged at step: 137, final residual: 1.72036e-08  Nonlinear convergence: ||u - u_old|| = 0.00314358
Linear solver converged at step: 250, final residual: 3.37166e-12  Nonlinear convergence: ||u - u_old|| = 1.96889e-05
 Nonlinear solver converged at step 3


*** Solving time step 9, time = 0.1 ***
Linear solver converged at step: 11, final residual: 0.00022001  Nonlinear convergence: ||u - u_old|| = 0.0505868
Linear solver converged at step: 69, final residual: 1.08054e-05  Nonlinear convergence: ||u - u_old|| = 0.0182682
Linear solver converged at step: 250, final residual: 3.49658e-08  Nonlinear convergence: ||u - u_old|| = 0.0127024
Linear solver converged at step: 250, final residual: 2.95441e-10  Nonlinear convergence: ||u - u_old|| = 5.53751e-05
 Nonlinear solver converged at step 3


*** Solving time step 10, time = 0.11 ***
Linear solver converged at step: 11, final residual: 0.000139574  Nonlinear convergence: ||u - u_old|| = 0.0322373
Linear solver converged at step: 104, final residual: 4.33263e-06  Nonlinear convergence: ||u - u_old|| = 0.0126853
Linear solver converged at step: 250, final residual: 4.09989e-08  Nonlinear convergence: ||u - u_old|| = 0.00565014
Linear solver converged at step: 250, final residual: 1.93796e-10  Nonlinear convergence: ||u - u_old|| = 5.64589e-05
 Nonlinear solver converged at step 3


*** Solving time step 11, time = 0.12 ***
Linear solver converged at step: 10, final residual: 0.000145884  Nonlinear convergence: ||u - u_old|| = 0.0205125
Linear solver converged at step: 48, final residual: 4.75062e-06  Nonlinear convergence: ||u - u_old|| = 0.00781029
Linear solver converged at step: 250, final residual: 1.24408e-08  Nonlinear convergence: ||u - u_old|| = 0.00689828
Linear solver converged at step: 250, final residual: 1.17963e-10  Nonlinear convergence: ||u - u_old|| = 2.01527e-05
 Nonlinear solver converged at step 3


*** Solving time step