Example problem: Adaptive solution of the 2D advection diffusion equation

In this example we discuss the adaptive solution of the 2D advection-diffusion problem

Two-dimensional advection-diffusion problem in a rectangular domain Solve
in the rectangular domain , with Dirichlet boundary conditions
where the |

We choose the forcing function and the boundary conditions such that

is the exact solution of the problem. For large values of the exact solution approaches a step, oriented at an angle against the axis.

In the computations we will impose the "wind"

illustrated in this vector plot:

The graph below shows a plot of the solution, computed at various levels of mesh adaptation, for and a Peclet number of

More interesting is the following plot which shows the solution for the same parameter values and boundary conditions, but for a zero forcing function,

The plot nicely illustrates the physical effects represented by the (unforced) advection diffusion equation. If represents the concentration of a chemical that is advected by the velocity field , while being dispersed by molecular diffusion, the advection-diffusion equation describes the steady-state concentration of this chemical. In this context the Peclet number is a measure of the relative importance of advective and diffusive effects. For very small Peclet number, the concentration is determined predominantly by diffusive effects – as , the advection diffusion equation approaches the Poisson equation. Conversely, at large values of the Peclet number, the concentration is determined predominantly by advective effects. The chemical is "swept along" by the flow and diffusive effects are only important in thin "boundary" or "shear" layers in which the concentration varies over short lengthscales. These can be seen clearly in the most finely resolved solution above.

Overall, the structure of the driver code is very similar to that in the corresponding Poisson example. The only difference is that we have to specify function pointers to the source and the "wind" functions, which are passed to the problem constructor. We create the problem, perform a self-test, set the global parameters that affect the solution and solve the problem using `oomph-lib's`

"black-box" adaptive Newton solver.

//===== start_of_main=====================================================

/// Driver code for 2D AdvectionDiffusion problem

//========================================================================

int main()

{

//Set up the problem

//------------------

// Create the problem with 2D nine-node refineable elements from the

// RefineableQuadAdvectionDiffusionElement family. Pass pointer to

// source and wind function.

RefineableAdvectionDiffusionProblem<RefineableQAdvectionDiffusionElement<2,

3> > problem(&TanhSolnForAdvectionDiffusion::source_function,

// Check if we're ready to go:

//----------------------------

cout << "\n\n\nProblem self-test ";

if (problem.self_test()==0)

{

cout << "passed: Problem can be solved." << std::endl;

}

else

{

throw OomphLibError("Self test failed",

OOMPH_CURRENT_FUNCTION,

OOMPH_EXCEPTION_LOCATION);

}

// Set the orientation of the "step" to 45 degrees

// Choose a large value for the steepness of the "step"

// Solve the problem, performing up to 4 adptive refinements

problem.newton_solve(4);

//Output the solution

problem.doc_solution();

} // end of main

The specification of the source function and the exact solution in the namespace `TanhSolnForAdvectionDiffusion`

is similar to that for the Poisson examples. The only difference is the inclusion of the Peclet number and the "wind" function.

//======start_of_namespace============================================

/// Namespace for exact solution for AdvectionDiffusion equation

/// with "sharp" step

//====================================================================

namespace TanhSolnForAdvectionDiffusion

{

/// Peclet number

double Peclet=200.0;

/// Parameter for steepness of step

double Alpha;

/// Parameter for angle of step

double TanPhi;

/// Exact solution as a Vector

{

}

/// Exact solution as a scalar

{

}

/// Source function required to make the solution above an exact solution

{

double x=x_vect[0];

double y=x_vect[1];

source =

}

/// Wind

{

wind[0]=sin(6.0*x[1]);

wind[1]=cos(6.0*x[0]);

}

} // end of namespace

The problem class is very similar to those used in the corresponding Poisson examples. The only change is that we use the function `Problem::actions_before_adapt()`

to document the progress of the automatic spatial adaptation. For this purpose, we store a `DocInfo`

as private member data in the `Problem`

. This allows us to increment the counter that labels the output files, accessible from `DocInfo::number()`

, whenever a new solution has been documented.

//====== start_of_problem_class=======================================

/// 2D AdvectionDiffusion problem on rectangular domain, discretised

/// with refineable 2D QAdvectionDiffusion elements. The specific type

/// of element is specified via the template parameter.

//====================================================================

template<class ELEMENT>

{

public:

/// Constructor: Pass pointer to source and wind functions

AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,

AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt);

/// Destructor. Empty

/// \short Update the problem specs before solve: Reset boundary conditions

/// to the values from the tanh solution.

void actions_before_newton_solve();

/// Update the problem after solve (empty)

void actions_after_newton_solve(){}

/// Actions before adapt: Document the solution

void actions_before_adapt()

{

// Doc the solution

doc_solution();

// Increment label for output files

Doc_info.number()++;

}

/// \short Doc the solution.

void doc_solution();

/// \short Overloaded version of the problem's access function to

/// the mesh. Recasts the pointer to the base Mesh object to

/// the actual mesh type.

RefineableRectangularQuadMesh<ELEMENT>* mesh_pt()

{

return dynamic_cast<RefineableRectangularQuadMesh<ELEMENT>*>(

Problem::mesh_pt());

}

private:

/// DocInfo object

DocInfo Doc_info;

/// Pointer to source function

AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt Source_fct_pt;

/// Pointer to wind function

AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt Wind_fct_pt;

}; // end of problem class

The constructor is practically identical to the constructors used in the various Poisson examples. We specify the output directory in the `Problem's`

`DocInfo`

object, create the mesh and an error estimator, and apply the boundary conditions by pinning the nodal values on the Dirichlet boundaries.

//=====start_of_constructor===============================================

/// \short Constructor for AdvectionDiffusion problem: Pass pointer to

/// source function.

//========================================================================

template<class ELEMENT>

AdvectionDiffusionEquations<2>::AdvectionDiffusionSourceFctPt source_fct_pt,

AdvectionDiffusionEquations<2>::AdvectionDiffusionWindFctPt wind_fct_pt)

: Source_fct_pt(source_fct_pt), Wind_fct_pt(wind_fct_pt)

{

// Set output directory

Doc_info.set_directory("RESLT");

// Setup mesh

// # of elements in x-direction

unsigned n_x=4;

// # of elements in y-direction

unsigned n_y=4;

// Domain length in x-direction

double l_x=1.0;

// Domain length in y-direction

double l_y=2.0;

// Build and assign mesh

Problem::mesh_pt() =

new RefineableRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);

// Create/set error estimator

mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;

// Set the boundary conditions for this problem: All nodes are

// free by default -- only need to pin the ones that have Dirichlet

// conditions here

unsigned num_bound = mesh_pt()->nboundary();

for(unsigned ibound=0;ibound<num_bound;ibound++)

{

unsigned num_nod= mesh_pt()->nboundary_node(ibound);

for (unsigned inod=0;inod<num_nod;inod++)

{

mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);

}

} // end loop over boundaries

We complete the problem setup by passing the function pointers to the source and wind functions, and the pointer to the Peclet number to the elements. Finally, we set up the equation numbering scheme.

// Complete the build of all elements so they are fully functional

// Loop over the elements to set up element-specific

// things that cannot be handled by the (argument-free!) ELEMENT

// constructor: Pass pointer to source function

unsigned n_element = mesh_pt()->nelement();

for(unsigned i=0;i<n_element;i++)

{

// Upcast from GeneralsedElement to the present element

ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));

//Set the source function pointer

el_pt->source_fct_pt() = Source_fct_pt;

//Set the wind function pointer

el_pt->wind_fct_pt() = Wind_fct_pt;

// Set the Peclet number

el_pt->pe_pt() = &TanhSolnForAdvectionDiffusion::Peclet;

}

// Setup equation numbering scheme

cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;

} // end of constructor

As before, we use the `Problem::actions_before_newton_solve()`

function to set/update the boundary conditions.

//========================================start_of_actions_before_newton_solve===

/// Update the problem specs before solve: (Re-)set boundary conditions

/// to the values from the tanh solution.

//========================================================================

template<class ELEMENT>

{

// How many boundaries are there?

unsigned num_bound = mesh_pt()->nboundary();

//Loop over the boundaries

for(unsigned ibound=0;ibound<num_bound;ibound++)

{

// How many nodes are there on this boundary?

unsigned num_nod=mesh_pt()->nboundary_node(ibound);

// Loop over the nodes on boundary

for (unsigned inod=0;inod<num_nod;inod++)

{

// Get pointer to node

Node* nod_pt=mesh_pt()->boundary_node_pt(ibound,inod);

// Extract nodal coordinates from node:

Vector<double> x(2);

x[0]=nod_pt->x(0);

x[1]=nod_pt->x(1);

// Compute the value of the exact solution at the nodal point

Vector<double> u(1);

// Assign the value to the one (and only) nodal value at this node

nod_pt->set_value(0,u[0]);

}

}

} // end of actions before solve

The function `doc_solution`

(...) is identical to that in the Poisson example. We output the solution, the exact solution and the error.

//===============start_of_doc=============================================

/// Doc the solution

//========================================================================

template<class ELEMENT>

{

ofstream some_file;

char filename[100];

// Number of plot points: npts x npts

unsigned npts=5;

// Output solution

//-----------------

sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),

Doc_info.number());

some_file.open(filename);

mesh_pt()->output(some_file,npts);

some_file.close();

// Output exact solution

//----------------------

sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),

Doc_info.number());

some_file.open(filename);

mesh_pt()->output_fct(some_file,npts,TanhSolnForAdvectionDiffusion::get_exact_u);

some_file.close();

// Doc error and return of the square of the L2 error

//---------------------------------------------------

double error,norm;

sprintf(filename,"%s/error%i.dat",Doc_info.directory().c_str(),

Doc_info.number());

some_file.open(filename);

mesh_pt()->compute_error(some_file,TanhSolnForAdvectionDiffusion::get_exact_u,

error,norm);

some_file.close();

// Doc L2 error and norm of solution

cout << "\nNorm of error : " << sqrt(error) << std::endl;

cout << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;

} // end of doc

- Explore the change in the character of the solution of the unforced problem when the Peclet number is slowly increased from 0 to 200, say. Note how at small Peclet number, strong diffusive effects smooth out the rapid spatial variations imposed by the boundary conditions. Conversely, at large values of the Peclet number, the behaviour is dominated by advective effects. As a result, in regions where the "wind" is directed into the domain, the value of set by the Dirichlet boundary conditions is "swept" into the domain. In regions where the "wind" is directed out of the domain, the value of "swept along" by the flow in the interior "clashes" with the value prescribed by the boundary conditions and the solution adjusts itself over a very short length scale, leading to the development of thin "boundary layers".
- Explore the character of the solution on coarse meshes at large and small Peclet numbers. Note how at large Peclet numbers the solution on the coarse meshes displays strong "wiggles" throughout the domain. These only disappear once the mesh adaptation fully resolves the regions of rapid variation. We will explore this issue further in another example.

- The source files for this tutorial are located in the directory:
demo_drivers/advection_diffusion/two_d_adv_diff_adapt/ - The driver code is:
demo_drivers/advection_diffusion/two_d_adv_diff_adapt/two_d_adv_diff_adapt.cc

A pdf version of this document is available.