Getting started

In this chapter, we explain how to solve the Poisson equation in 1d, 2d and 3d. We also introduce every concept/object used.

Let \Omega \subset \mathbb{R}^d be a computational domain that is the image of a logical domain \mathcal{P}, i.e. a unit line (in 1d), square (in 2d) or a cube (in 3d) with a mapping function

\Omega = \mathcal{G} (\mathcal{P})

We consider the Poisson problem with Homogeneous Dirichlet boundary conditions:

\mbox{Find} ~  u \in H^1_0(\Omega), ~ \mbox{such that}
\\
- \nabla^2 u = f, ~~ \Omega

Using a Galerkin-Rietz (or more generally Galerkin-Petrov) method, we introduce:

  1. a discrete finite elements space \mathcal{V}_h = \mathbf{span}\{ \phi_j, 1 \leq j \leq n_V \} for trial functions
  2. a discrete finite elements space \mathcal{W}_h = \mathbf{span}\{ \psi_i, 1 \leq i \leq n_W \} for test functions (in the case of Galerkin-Rietz method, we have \mathcal{W}_h = \mathcal{V}_h )

We will skip details about the construction of such spaces. Please take a look at DISCO documentation for more details.

Let u_h \in \mathcal{V}_h such that u_h = \sum_{j=1}^{n_V} u_j \phi_j. Then the weak formulation for the Poisson equation writes

\sum_{j=1}^{n_V} u_j \int_{\Omega} \nabla \phi_j \cdot \nabla \psi_i = \int_{\Omega} f \psi_i, \quad \forall  1 \leq i \leq n_W

Now, because we want to use the Galerkin-Rietz approach, we have

\sum_{j=1}^{n} u_j \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i = \int_{\Omega} f \phi_i, \quad \forall  1 \leq i \leq n

which can be written in a matrix form

M U = F

where U denotes the coefficients (u_j, ~ 1 \leq j \leq n) and F is a vector consisting of the terms \int_{\Omega} f \phi_i for 1 \leq i \leq n. Finally, the matrix M is called the stiffness matrix and its entries are defined as

M_{ij} = \int_{\Omega} \nabla \phi_j \cdot \nabla \phi_i

For the sake of consistency with FEMA we will denote our basis function by v_i and v_j rather than \phi_i and \phi_j. In this case, in order to solve the Poisson problem, one needs to

  1. Assemble the matrix M of entries

M_{ij} = \int_{\Omega} \nabla v_j \cdot \nabla v_i

  1. Assemble the right hand side F of entries

F_{i} = \int_{\Omega} \nabla f v_i

  1. Solve the linear system

M U = F

The aim of this tutorial, is to show you how to solve numerically this problem, using the different concepts and libraries provided within the JOREK-Django framework.

Poisson 1d

In this section, we’ll be considering the file ex01_01.F90 that you can find in fema/tutorials/examples/1d/poisson/.

How to compile the example

As explained before, you need to create a build directory:

> cd
> mkdir build_poisson_1d && cd build_poisson_1d

Configuring and compiling the simulation is done by calling the make_project.sh script:

> $CLAPP_DIR/scripts/make_project.sh PATH_TO_FEMA/tutorials/examples/1d/poisson/

Note

Notice that we are in fact compiling other files too.

How to run the example

First, copy the following input directory:

> cp -R $CLAPP_DIR/inputs/fema/1d/poisson/dirichlet_ex1 input

You can then run the example:

> bin/poisson_1d_ex01_01 input/parameters_sim_1d.nml input/solver_driver.nml

Before going forward, let us take a short break and see what is inside these files. If you run:

> ls input

You’ll see the following files:

parameters_bspline.nml
parameters_ddm.nml
parameters_sim_1d.nml
solver_driver.nml
solver_lapack_lu.nml

So why so many files, while we are only calling the binary with two files? let’s first open the file parameters_sim_1d.nml:

> more parameters_sim_1d.nml
&ASSEMBLY
  ID=        1,
 /
&DDM
  DDM_NAME="default",
  DDM_FILENAME="input/parameters_ddm.nml",
/
&DISCRETIZATION
 NAME_1="bspline",
 FILENAME_1="input/parameters_bspline.nml",
 /

There are 3 sections

  1. ASSEMBLY: we specify here which assembler object we want to use for our weak formulations. In the 1d case, there’s only one assembler object, and we refeer to it by the ID 1.
  2. DDM: this sections concerns the domain decomposition method we’re using. There are two entries:
  • DDM_NAME: specifies the type of the ddm object to allocate
  • DDM_FILENAME: the filename describing the parameters for the ddm.

Basically, the first attribut allows us to know which type to allocate then we create it from file, using the read_from_file method.

  1. DISCRETIZATIONS: this sections concerns the kind of discretization that will be used in our simulation. Again, there are two entries:
  • NAME_1: specifies the type of the finite elements space to allocate
  • FILENAME_1: specifies the input file describing the descritization, following the same idea as for the ddm.

So, in this example, we are using a B-Spline discretization, with a default ddm.Next, let’s see what’s inside the two files:

> more input/parameters_ddm.nml
&P_DIMENSION
  DDM_DIM          = 1,
 /
&PARAMETERS
  PARTITIONER_NAME = 'DEFAULT',
  N_PROCS          = 1,
 /

So if, you want to run the simulation in parallel, you’ll only need to change the value of N_PROCS to match the number of processors you’re using when invoking MPI.

Now, let’s take a look to the other file:

> more input/parameters_bspline.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 64,
 DEGREES=          5,
 N_PROCS=          1,
 BC_MIN=           0,
 BC_MAX=           0,
 /

So as you can see, we are using a B-Spline finite elements discrezation on a grid of n_elements = 64 with a polynomial degree 5.

Finally, the boundary conditions are specified by giving the code ID 0 which means that we are removing the interpolatory basis functions from the finite elements space.

Note

As you can see, there is some redundency; for instance, we are specifying that the dimension of our problem is 1 in two files. The reason is to preserve consistency for our objects.

Note

N_PROCS is also an attribute for the discretization parameters file. In fact, if you use the static type related to the B-Spline discretization, then you don’t need to have a simulation or a ddm parameter files.

The last argument of our binary file is input/solver_driver.nml. Let’s see what it contains:

> more input/solver_driver.nml
&linear_solver

SOLVER_NAME="lapack_lu"
SOLVER_FILENAME="input/solver_lapack_lu.nml",

/

As for the simulation parameter file, we see that we have two arguments. We first give the name of the solver SOLVER_NAME needed to allocate the corresponding type. Then we give the file that contains the parameters for the solver SOLVER_FILENAME. Here for instance:

> more input/solver_lapack_lu.nml
&linear_solver

INT_VERBOSE=1
INT_TRANSPOSE=0

/

Nothing special as we are using a direct solver from LAPACK.

Note

for more details about the available linear solvers, we refer the reader to PLAF documentation.

Note

the driver linear solver allows the user to specify the linear solver to use at the run time.

Now let’s run our binary:

> bin/poisson_1d_ex01_01 input/parameters_sim_1d.nml input/solver_driver.nml
 >>>> context
 * assembly id            1
 * dimension            1
 * ddm name default
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 <<<<
 >>>> context
 * assembly id            1
 * dimension            1
 * ddm name default
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 <<<<
 >>>>> parameters
 partitioner_id   :            0
 partitioner_name : DEFAULT
 ddm_dim          :            1
 n_procs         :            1
 <<<<<
Assembly time =  0.000 seconds, on         64 elements.
Assembly time =  0.004 seconds, on         64 elements.
 >> norms :
   4.2632886523999863E-007   1.3249946937903078E-005

Visualization

We provide two ways for the visualization. The first one is based on the output object offred by FEMA which will be detailed in the next part. The result of the output export is a tecplot file that can be read by a visualization tool like VisIt or ParaView.

The second one is based on plot_field script available in DISCO. In the following example, we’ll directly export the field and the mapping into namelist files that we’ll process using Python. After executing our binary, we see the following new files:

phi.nml
mapping.nml

Now in order to visualize the solution stored in the file phi.nml, run:

> python $CLAPP_DIR/scripts/disco/plot_field.py phi.nml mapping.nml

The result is the following plot

alternate text

Numerical solution u_h

Dive into the code

It’s time to open the ex01_01.F90 file.

Initialize and Finalize subroutines

Before going further, let’s take a look to the first and last calls within our program:

! ...
call jrk_initialize(i_err)
! ...

and

! ...
call jrk_finalize(i_err)
! ...

These subroutines are used to initialize and finalize the simulation. For instance, they include calls to MPI if used.

Declarations and polymorphic variables

Let’s take a look to our declarations. First you notice that the objects trial_space and test_space are declared as abstract

class(dsc_t_space_abstract), allocatable :: trial_space
class(dsc_t_space_abstract), allocatable :: test_space

This is also the case for the field phi

class(dsc_t_field_abstract), allocatable :: phi

and the assembler object

class(jrk_t_assembler_abstract),  allocatable :: assembler

In fact, we allocate the type of these objects dynamically, using

! ... allocates and creates the finite elements space from the context
call jrk_allocate_space (trial_space, trial_context)
call jrk_allocate_space (test_space, test_context)
! ...

and

! ... allocates the field type from the context
call jrk_allocate_field (phi, trial_context)
! ...

and

! ... allocates the assembler type from the context
call jrk_allocate_assembler (assembler, trial_context)
! ...

Todo

maybe it should be better to have one single subroutine for allocation jrk_allocate that will call jrk_allocate_space, jrk_allocate_field etc

Please notice that the subroutine jrk_allocate_space is doing two things. First it allocates the right type for the space, then it creates it. On the other hand, invoking jrk_allocate_field or jrk_allocate_assembler are only doing the allocation. Therefor, the user must call the create method for the field phi

! ... creates field on the trial space
call phi % create(space=trial_space)
! ...

Then for the assembler

! ... creates an assembler object for the couple (test_space, trial_space) with
!       1 field
!       1 rhs
!       1 matrix
!       a vector to store the norms
!       a mapping (otherwise, the mapping associated to the trial space is used)
!       the ddm parameters
call assembler % create( trial_space                 &
                     & , test_space                  &
                     & , n_fields   = 1              &
                     & , n_rhs      = 1              &
                     & , n_matrices = 1              &
                     & , norms = norms               &
                     & , mapping = mapping           &
                     & , ddm_parameters = ddm_params )
! ...

Now, let’s go back to our declarations

class(dsc_t_space_abstract),      allocatable :: trial_space
class(dsc_t_space_abstract),      allocatable :: test_space
class(dsc_t_field_abstract),      allocatable :: phi
class(jrk_t_assembler_abstract),  allocatable :: assembler
type(jrk_t_context),              target      :: trial_context
type(jrk_t_context),              target      :: test_context
type(plf_t_vector),               target      :: rhs
type(plf_t_matrix_csr),           target      :: matrix
type(plf_t_vector),               target      :: norms
type(spl_t_mapping_1d),           target      :: mapping
type(plf_t_ddm_parameters),       target      :: ddm_params
type(plf_t_linear_solver_driver), target      :: linear_solver

You see that the only polymorphic objects are the one we already introduced: trial_space, test_space and phi. The remaining objects are defined using the keyword type and not class.

Geometry and Mapping

In this example we are using a linear mapping between the default logical patch (a unit line [0,1] ) and the line [P_1,P_2]

! ... creates a linear mapping between (0,1) and (P_1, P_2)
P_1 = (/ -2.0_plf_rk /)
P_2 = (/  2.0_plf_rk /)

call spl_mapping_linear(mapping, P_1, P_2)
call mapping % export('mapping.nml')
! ...

Maybe you noticed that we created our finite elements space without specifying the mapping. There are in fact three ways to use the mapping:

  1. create the space with the mapping.
  2. create the space then associate the mapping using the method set_mapping.
  3. create the space without a mapping, but call the assembler with the argument mapping.

In this example we choosed the 3 approach.

DDM

The ddm concept allows the partitioning of our discrete space, matrices, ... As shown before, the ddm parameter file is given inside the context one. We can then create a ddm parameter object using

! ... reads the ddm parameters
call ddm_params % read_from_file(trial_context % ptr_parameters % ddm_filename)
call ddm_params % print_info()
! ...

Linear Algebra

In order to solve the Poisson problem, we need to declare

  1. one matrix associated to one variable
  2. one rhs associated to one variable

This is done by the following lines

! ... allocation is postponed until we give the ddm
call matrix % create(n_block_rows=1, n_block_cols=1)
! ...

! ... allocation is postponed until we give the ddm
call rhs % create(n_blocks=1)
! ...

The number of blocks (rows and cols) is in fact to the number of variables per degree of freedom.

More about the assembler object

Now, we have our matrix and rhs object, but how to link them to the assembler object?

The assembler offers three methods to specify and associate objects like

  • matrices using the method set_matrix
  • rhs using the method set_rhs
  • fields using the method set_field

Please notice that these methods take an optional id that describes its position within the assembler data structure.

! ... associates the matrix, rhs and field to the assembler
call assembler % set_matrix(matrix, i_matrix=1)
call assembler % set_rhs(rhs, i_rhs=1)
call assembler % set_field(phi, i_field=1)
! ...

Now that we have associated all these objects, the only thing that remains is to specify the weak formulation and then assemble it.

Because we want to assemble two objects (matrix and rhs) we need to associate the assembler to two subroutines

  • one for the matrix and computes the integral between the trial and test basis functions v_j and v_i. This is done inside the subroutine build_matrix_stiffness.
  • one for the rhs and computes the integral between the function f and test basis function v_i. This is done inside the subroutine build_rhs.
! ... sets the weak formulation for both the matrix and rhs
assembler % ptr_matrix_contribution => build_matrix_stiffness
assembler % ptr_rhs_contribution    => build_rhs

call assembler % assemble(mapping=mapping)
! ...

Now let’s take a look to these subroutines. First, you’ll notice that build_matrix_stiffness is not defined in our example. In fact, FEMA provides some usefull weak formulations in the gallery package (mass, stiffness for example). You can find this subroutine in the file fema/fortran/src/gallery/gallery_1d.F90

! ..................................................
subroutine build_matrix_stiffness(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(jrk_rk) :: vi_0
  real(jrk_rk) :: vi_x
  real(jrk_rk) :: vj_0
  real(jrk_rk) :: vj_x
  real(jrk_rk) :: wvol

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % matrix_contribution = 0.0_jrk_rk
  ! ...

  ! ...
  do i_point = 1, n_points

    vi_0  = self % arr_vi(1, 1, i_point)
    vi_x  = self % arr_vi(1, 2, i_point)

    vj_0  = self % arr_vj(1, 1, i_point)
    vj_x  = self % arr_vj(1, 2, i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % matrix_contribution = self % matrix_contribution + vi_x * vj_x * wvol
  end do
  ! ...

end subroutine build_matrix_stiffness
! ..................................................

You recover here our naming convention on the basis functions. You also notice how we compute the derivatives (vi_x for example). The aim of this subroutine (which is called for every couple (v_j, v_i) is to evaluate and cumulate the values of basis functions, their derivatives and some input data on the quadrature points, then update the attribut matrix_contribution (more details about this attribut will be given later).

Now, for the rhs, it is defined in our poisson file

! ..................................................
subroutine build_rhs(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: f_value
  real(plf_rk) :: vi_0
  real(plf_rk) :: vi_x
  real(plf_rk) :: wvol
  real(plf_rk) :: x
  real(kind=plf_rk), dimension(self % n_points) :: f_values

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % rhs_contribution = 0.0_plf_rk
  ! ...

  ! ...
  f_values = (k1**2) * sin(k1 * self % arr_x (1, :))
  ! ...

  ! ...
  do i_point = 1, n_points

    x    = self % arr_x (1, i_point)

    vi_0 = self % arr_vi(1, 1, i_point)
    vi_x = self % arr_vi(1, 2, i_point)

    f_value = f_values(i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % rhs_contribution = self % rhs_contribution + f_value * vi_0 * wvol
  end do
  ! ...

end subroutine build_rhs
! ..................................................

We recover the same accumulation principal that was introduced for the matrix. Here, the aim is to update the attribut rhs_contribution.

Linear Algebra Cont’

At this point, we have assembled our matrix and rhs. Let’s now solve the associated linear system. For this purpose, we need to create our linear solver (remember it was a driver one)

! ... creates a driver linear solver
call linear_solver % create(matrix, filename_solver)
! ...

Solving the linear system is done by the following lines

! ... solves and sets the field
allocate(x(trial_space % ptr_numbering % n_vertex))
allocate(y(test_space % ptr_numbering % n_vertex))

x = 0.0_plf_rk
y = 0.0_plf_rk
call rhs % get(y)
call linear_solver % solve(y, x)
call phi % set(x)
! ...

Computing the L^2 and H^1 norms

Following the same idea for the weak formulation, we will associate to the assembler object the following two subroutines

! ... associates fields and norms evaluation
assembler % ptr_fields_evaluation   => evaluate_fields
assembler % ptr_norms_evaluation    => evaluate_norms

call assembler % assemble(mapping=mapping)
! ...

The first subroutine evaluate_fields is used here to evaluate u_h (its first order derivative) on the quadrature points.

Note

This subroutine can be used if one deals with a nonlinear problem.

evaluate_fields is defined as follows

! ..................................................
subroutine evaluate_fields(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_trial_basis
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: vj_0
  real(plf_rk) :: vj_x

  ! ...
  do i_point = 1, self % n_points

    vj_0  = self % arr_vj(1, 1, i_point)
    vj_x  = self % arr_vj(1, 2, i_point)

    self % field_values(1, 1, :, i_point) = self % field_values(1, 1, :, i_point) &
                                  & + vj_0 * self % field_coeffs(:)

    self % field_values(1, 2, :, i_point) = self % field_values(1, 2, :, i_point) &
                                  & + vj_x * self % field_coeffs(:)
  end do
  ! ...

end subroutine evaluate_fields
! ..................................................

Lets’ take a look at the definition of field_values. It is allocated in assembler_1d.F90 as follows

allocate(self % field_values (self % c_dim, 2, self % n_fields, self % n_points))

Therefor,

  • first index for the dimension of the field (here it’s a scalar, then it’s 1)
  • second index for the derivative order of the field (here it’s 2)
  • third index for the field id (here it’s 1)
  • fourth index describes the quadrature points.

Finally the subroutine evaluate_norms is defined as

! ..................................................
subroutine evaluate_norms(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_diagnostic
  integer :: i_field
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: wvol
  real(plf_rk) :: norm_l2
  real(plf_rk) :: norm_h1
  real(plf_rk) :: f_0
  real(plf_rk) :: f_x
  real(plf_rk) :: field_0
  real(plf_rk) :: field_x
  real(kind=plf_rk), dimension(2, self % n_points) :: f_values

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  f_values(1, :) = sin(k1 * self % arr_x (1, :))
  f_values(2, :) = k1 * cos(k1 * self % arr_x (1, :))
  ! ...

  ! ...
  i_diagnostic = 0
  do i_field = 1, self % n_fields
    ! ...
    norm_l2 = 0.0_plf_rk
    norm_h1 = 0.0_plf_rk
    do i_point = 1, n_points

      f_0 = f_values(1, i_point)
      f_x = f_values(2, i_point)

      field_0 = self % field_values(1, 1, i_field, i_point)
      field_x = self % field_values(1, 2, i_field, i_point)

      wvol = self % weights(i_point) * self % jacobians(i_point)

      norm_l2 = norm_l2 + (field_0 - f_0)**2 * wvol
      norm_h1 = norm_h1 + (field_x - f_x)**2 * wvol
    end do
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_l2)
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_h1)
    ! ...
  end do
  ! ...

end subroutine evaluate_norms
! ..................................................

You notice that we are using the field_values array and then we update the pointer self % ptr_norms to the norms vector that we gave when creating the assembler object.

Output plot

The output object offred by FEMA allows to evaluate a field on the corresponding mesh, providing a given number of points per element.

! ..................................................
call output % create( export_name = "poisson_1d_ex01_01", &
                      & space     = trial_space , &
                      & n_fields  = 1           ,     &
                      & n_points_per_elements = 2 )


call output % set_field(phi, i_field=1)
call output % export()
! ..................................................

In this example, the export will generate the file poisson_1d_ex01_01_0.tp.

Memory deallocation

Don’t forget that we need to destroy and deallocate the memory at the end of our simulation

! ... deallocates memory
call linear_solver % free()
call assembler % free()
call output % free()
call norms % free()
call rhs % free()
call phi % free()
call test_space % free()
call trial_space % free()
call mapping % free()
call trial_context % free()
call test_context % free()

deallocate(trial_space)
deallocate(test_space)
deallocate(phi)
deallocate(assembler)
! ...

Poisson 2d

In this section, we’ll be considering the ex01_01.F90 that you can find in fema/tutorials/examples/2d/poisson/.

How to compile the example

As explained before, you need to create a build directory:

> cd
> mkdir build_poisson_2d && cd build_poisson_2d

Configuring and compiling the simulation is done by calling the make_project.sh script:

> $CLAPP_DIR/scripts/make_project.sh PATH_TO_FEMA/tutorials/examples/2d/poisson/

Note

Notice that we are in fact compiling other files too.

How to run the example

First, copy the following input directory:

> cp -R $CLAPP_DIR/inputs/fema/2d/poisson/dirichlet_ex1 input

You can then run the example:

> bin/poisson_2d_ex01_01 input/parameters_sim_1d1d.nml input/solver_driver.nml

Before going forward, let us take a short break and see what is inside these files. If you run:

> ls input

You’ll see the following files:

parameters_bspline_u.nml
parameters_bspline_v.nml
parameters_ddm.nml
parameters_sim_1d1d.nml
solver_driver.nml
solver_lapack_lu.nml

So why so many files, while we are only calling the binary with two files? let’s first open the file parameters_sim_1d1d.nml:

> more parameters_sim_1d1d.nml
&ASSEMBLY
  ID=        11,
 /
&DDM
  DDM_NAME="tensor",
  DDM_FILENAME="input/parameters_ddm.nml",
/
&DISCRETIZATION
  NAME_1="bspline",
  FILENAME_1="input/parameters_bspline_u.nml",
  NAME_2="bspline",
  FILENAME_2="input/parameters_bspline_v.nml",
 /

There are 3 sections

  1. ASSEMBLY: we specify here which assembler object we want to use for our weak formulations. In the 2d case, there are two assembler objects (the 2d and 1d x 1d tensor). For the tensor case, the ID 11.
  2. DDM: this sections concerns the domain decomposition method we’re using. There are two entries:
  • DDM_NAME: specifies the type of the ddm object to allocate
  • DDM_FILENAME: the filename describing the parameters for the ddm.

Basically, the first attribut allows us to know which type to allocate then we create it from file, using the read_from_file method.

  1. DISCRETIZATIONS: this sections concerns the kind of discretization that will be used in our simulation. Again, there are four entries:
  • NAME_1: specifies the type of the finite elements space to allocate in the u-direction
  • FILENAME_1: specifies the input file describing the descritization, following the same idea as for the ddm in the u-direction.
  • NAME_2: specifies the type of the finite elements space to allocate in the v-direction
  • FILENAME_2: specifies the input file describing the descritization, following the same idea as for the ddm in the v-direction.

So, in this example, we are using a B-Spline discretization, with a tensor ddm.Next, let’s see what’s inside the two files:

> more input/parameters_ddm.nml
&P_DIMENSION
  DDM_DIM          = 2,
 /
&PARAMETERS
  PARTITIONER_NAME = 'TENSOR',
  N_PROCS          = 1, 1,
 /

So if, you want to run the simulation in parallel, you’ll only need to change the values of N_PROCS to match the total number of processors you’re using when invoking MPI.

Now, let’s take a look to the other files:

> more input/parameters_bspline_u.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 32,
 DEGREES=          3,
 BC_MIN=           0,
 BC_MAX=           0,
 /

> more input/parameters_bspline_v.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 32,
 DEGREES=          3,
 BC_MIN=           0,
 BC_MAX=           0,
 /

So as you can see, we are using a B-Spline finite elements discrezation on a grid of n_elements = 32 \times 32 with a polynomial degrees (3,3).

Finally, the boundary conditions are specified by giving the code ID 0 which means that we are removing the interpolatory basis functions from the finite elements space.

Now let’s run our binary:

> bin/poisson_2d_ex01_01 input/parameters_sim_1d1d.nml input/solver_driver.nml
 >>>> context
 * assembly id           11
 * dimension            2
 * ddm name tensor
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 * discretization_2 bspline
 <<<<
 >>>> context
 * assembly id           11
 * dimension            2
 * ddm name tensor
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 * discretization_2 bspline
 <<<<
 >>>>> parameters
 partitioner_id   :            1
 partitioner_name : TENSOR
 ddm_dim          :            2
 n_procs         :            1           1
 <<<<<
Assembly time =  0.276 seconds, on       1024 elements.
Assembly time =  0.040 seconds, on       1024 elements.
 >> norms :
   1.2368090699991091E-003   5.7994854538192669E-002

Visualization

We provide two ways for the visualization. The first one is based on the output object offred by FEMA which will be detailed in the next part. The result of the output export is a tecplot file that can be read by a visualization tool like VisIt or ParaView.

The second one is based on plot_field script available in DISCO. In the following example, we’ll directly export the field and the mapping into namelist files that we’ll process using Python. After executing our binary, we see the following new files:

phi.nml
mapping.nml

Now in order to visualize the solution stored in the file phi.nml, run:

> python $CLAPP_DIR/scripts/disco/plot_field.py phi.nml mapping.nml

The result is the following plot

alternate text

Numerical solution u_h

Dive into the code

It’s time to open the ex01_01.F90 file.

Initialize and Finalize subroutines

This part is exactly similar to the 1d case.

Declarations and polymorphic variables

This part is exactly similar to the 1d case.

Geometry and Mapping

In this example we are using a bilinear mapping between the default logical patch (a unit square [0,1] \times [0,1] ) and the quadrangle [[P_{11},P_{12}], [P_{21},P_{22}]]

! ... creates a bilinear mapping
P_11 = (/ -2.0_jrk_rk, -2.0_jrk_rk /)
P_21 = (/  2.0_jrk_rk, -2.0_jrk_rk /)
P_12 = (/ -2.0_jrk_rk,  2.0_jrk_rk /)
P_22 = (/  2.0_jrk_rk,  2.0_jrk_rk /)

call spl_mapping_bilinear(mapping, P_11, P_12, P_21, P_22)
call mapping % export('mapping.nml')
! ...

Maybe you noticed that we created our finite elements space without specifying the mapping. There are in fact three ways to use the mapping:

  1. create the space with the mapping.
  2. create the space then associate the mapping using the method set_mapping.
  3. create the space without a mapping, but call the assembler with the argument mapping.

In this example we choosed the 3 approach.

DDM

This part is exactly similar to the 1d case.

Linear Algebra

This part is exactly similar to the 1d case.

More about the assembler object

Now, we have our matrix and rhs object, but how to link them to the assembler object?

The assembler offers three methods to specify and associate objects like

  • matrices using the method set_matrix
  • rhs using the method set_rhs
  • fields using the method set_field

Please notice that these methods take an optional id that describes its position within the assembler data structure.

! ... associates the matrix, rhs and field to the assembler
call assembler % set_matrix(matrix, i_matrix=1)
call assembler % set_rhs(rhs, i_rhs=1)
call assembler % set_field(phi, i_field=1)
! ...

Now that we have associated all these objects, the only thing that remains is to specify the weak formulation and then assemble it.

Because we want to assemble two objects (matrix and rhs) we need to associate the assembler to two subroutines

  • one for the matrix and computes the integral between the trial and test basis functions v_j and v_i. This is done inside the subroutine build_matrix_stiffness.
  • one for the rhs and computes the integral between the function f and test basis function v_i. This is done inside the subroutine build_rhs.
! ... sets the weak formulation for both the matrix and rhs
assembler % ptr_matrix_contribution => build_matrix_stiffness
assembler % ptr_rhs_contribution    => build_rhs

call assembler % assemble(mapping=mapping)
! ...

Now let’s take a look to these subroutines. First, you’ll notice that build_matrix_stiffness is not defined in our example. In fact, FEMA provides some usefull weak formulations in the gallery package (mass, stiffness for example). You can find this subroutine in the file fema/fortran/src/gallery/gallery_2d.F90

! ..................................................
subroutine build_matrix_stiffness(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(jrk_rk) :: vi_0
  real(jrk_rk) :: vi_x
  real(jrk_rk) :: vi_y
  real(jrk_rk) :: vj_0
  real(jrk_rk) :: vj_x
  real(jrk_rk) :: vj_y
  real(jrk_rk) :: wvol

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % matrix_contribution = 0.0_jrk_rk
  ! ...

  ! ...
  do i_point = 1, n_points

    vi_0  = self % arr_vi(1, 1, i_point)
    vi_x  = self % arr_vi(1, 2, i_point)
    vi_y  = self % arr_vi(1, 3, i_point)

    vj_0  = self % arr_vj(1, 1, i_point)
    vj_x  = self % arr_vj(1, 2, i_point)
    vj_y  = self % arr_vj(1, 3, i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % matrix_contribution = self % matrix_contribution &
                             & + ( vi_x * vj_x + vi_y * vj_y ) * wvol
  end do
  ! ...

end subroutine build_matrix_stiffness
! ..................................................

You recover here our naming convention on the basis functions. You also notice how we compute the derivatives (vi_x for example). The aim of this subroutine (which is called for every couple (v_j, v_i) is to evaluate and cumulate the values of basis functions, their derivatives and some input data on the quadrature points, then update the attribut matrix_contribution (more details about this attribut will be given later).

Now, for the rhs, it is defined in our poisson file

! ..................................................
subroutine build_rhs(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: f_value
  real(plf_rk) :: vi_0
  real(plf_rk) :: wvol
  real(kind=plf_rk), dimension(self % n_points) :: f_values

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % rhs_contribution = 0.0_plf_rk
  ! ...

  ! ...
  f_values = (k1**2 + k2**2) &
         & * sin(k1 * self % arr_x (1, :)) &
         & * sin(k2 * self % arr_x (2, :))
  ! ...

  ! ...
  do i_point = 1, n_points
    vi_0  = self % arr_vi(1, 1, i_point)

    f_value = f_values(i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % rhs_contribution = self % rhs_contribution + f_value * vi_0 * wvol
  end do
  ! ...

end subroutine build_rhs
! ..................................................

We recover the same accumulation principal that was introduced for the matrix. Here, the aim is to update the attribut rhs_contribution.

Linear Algebra Cont’

This part is exactly similar to the 1d case.

Computing the L^2 and H^1 norms

Following the same idea for the weak formulation, we will associate to the assembler object the following two subroutines

! ... associates fields and norms evaluation
assembler % ptr_fields_evaluation   => evaluate_fields
assembler % ptr_norms_evaluation    => evaluate_norms

call assembler % assemble(mapping=mapping)
! ...

The first subroutine evaluate_fields is used here to evaluate u_h (its first order derivative) on the quadrature points.

Note

This subroutine can be used if one deals with a nonlinear problem.

evaluate_fields is defined as follows

! ..................................................
subroutine evaluate_fields(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_trial_basis
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: vj_0
  real(plf_rk) :: vj_x
  real(plf_rk) :: vj_y

  ! ...
  do i_point = 1, self % n_points

    vj_0  = self % arr_vj(1, 1, i_point)
    vj_x  = self % arr_vj(1, 2, i_point)
    vj_y  = self % arr_vj(1, 3, i_point)

    self % field_values(1, 1, :, i_point) = self % field_values(1, 1, :, i_point) &
                                  & + vj_0 * self % field_coeffs(:)

    self % field_values(1, 2, :, i_point) = self % field_values(1, 2, :, i_point) &
                                  & + vj_x * self % field_coeffs(:)

    self % field_values(1, 3, :, i_point) = self % field_values(1, 3, :, i_point) &
                                  & + vj_y * self % field_coeffs(:)
  end do
  ! ...

end subroutine evaluate_fields
! ..................................................

Finally the subroutine evaluate_norms is defined as

! ..................................................
subroutine evaluate_norms(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_diagnostic
  integer :: i_field
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: wvol
  real(plf_rk) :: norm_l2
  real(plf_rk) :: norm_h1
  real(plf_rk) :: f_0
  real(plf_rk) :: f_x
  real(plf_rk) :: f_y
  real(plf_rk) :: field_0
  real(plf_rk) :: field_x
  real(plf_rk) :: field_y

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  i_diagnostic = 0
  do i_field = 1, self % n_fields
    ! ...
    norm_l2 = 0.0_plf_rk
    norm_h1 = 0.0_plf_rk
    do i_point = 1, n_points

      f_0 =      sin(k1 * self % arr_x (1, i_point)) * sin(k2 * self % arr_x (2, i_point))
      f_x = k1 * cos(k1 * self % arr_x (1, i_point)) * sin(k2 * self % arr_x (2, i_point))
      f_y = k2 * sin(k1 * self % arr_x (1, i_point)) * cos(k2 * self % arr_x (2, i_point))

      field_0 = self % field_values(1, 1, i_field, i_point)
      field_x = self % field_values(1, 2, i_field, i_point)
      field_y = self % field_values(1, 3, i_field, i_point)

      wvol = self % weights(i_point) * self % jacobians(i_point)

      norm_l2 = norm_l2 + (field_0 - f_0)**2 * wvol
      norm_h1 = norm_h1 + ((field_x - f_x)**2 + (field_y - f_y)**2) * wvol
    end do
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_l2)
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_h1)
    ! ...
  end do
  ! ...

end subroutine evaluate_norms
! ..................................................

You notice that we are using the field_values array and then we update the pointer self % ptr_norms to the norms vector that we gave when creating the assembler object.

Output plot

The output object offred by FEMA allows to evaluate a field on the corresponding mesh, providing a given number of points per element.

! ..................................................
call output % create( export_name = "poisson_2d_ex01_01", &
                      & space     = trial_space , &
                      & n_fields  = 1           ,     &
                      & n_points_per_elements = 2 )


call output % set_field(phi, i_field=1)
call output % export()
! ..................................................

In this example, the export will generate the file poisson_2d_ex01_01_0.tp.

Memory deallocation

This part is exactly similar to the 1d case.

Poisson 3d

In this section, we’ll be considering the ex01_01.F90 that you can find in fema/tutorials/examples/3d/poisson/.

How to compile the example

As explained before, you need to create a build directory:

> cd
> mkdir build_poisson_3d && cd build_poisson_3d

Configuring and compiling the simulation is done by calling the make_project.sh script:

> $CLAPP_DIR/scripts/make_project.sh PATH_TO_FEMA/tutorials/examples/3d/poisson/

Note

Notice that we are in fact compiling other files too.

How to run the example

First, copy the following input directory:

> cp -R $CLAPP_DIR/inputs/fema/3d/poisson/dirichlet_ex1 input

You can then run the example:

> bin/poisson_3d_ex01_01 input/parameters_sim_1d1d1d.nml input/solver_driver.nml

Before going forward, let us take a short break and see what is inside these files. If you run:

> ls input

You’ll see the following files:

parameters_bspline_u.nml
parameters_bspline_v.nml
parameters_bspline_w.nml
parameters_ddm.nml
parameters_sim_1d1d1d.nml
solver_driver.nml
solver_lapack_lu.nml

So why so many files, while we are only calling the binary with two files? let’s first open the file parameters_sim_1d1d1d.nml:

> more parameters_sim_1d1d1d.nml
&ASSEMBLY
  ID=        111,
 /
&DDM
  DDM_NAME="tensor",
  DDM_FILENAME="input/parameters_ddm.nml",
/
&DISCRETIZATION
  NAME_1="bspline",
  FILENAME_1="input/parameters_bspline_u.nml",
  NAME_2="bspline",
  FILENAME_2="input/parameters_bspline_v.nml",
  NAME_3="bspline",
  FILENAME_3="input/parameters_bspline_w.nml",
 /

There are 3 sections

  1. ASSEMBLY: we specify here which assembler object we want to use for our weak formulations. In the 3d case, there’s only one assembler object which corresponds to the 1d x 1d x 1d tensor case. We refeer to it by the ID 111.
  2. DDM: this sections concerns the domain decomposition method we’re using. There are two entries:
  • DDM_NAME: specifies the type of the ddm object to allocate
  • DDM_FILENAME: the filename describing the parameters for the ddm.

Basically, the first attribut allows us to know which type to allocate then we create it from file, using the read_from_file method.

  1. DISCRETIZATIONS: this sections concerns the kind of discretization that will be used in our simulation. Again, there are six entries:
  • NAME_1: specifies the type of the finite elements space to allocate in the u-direction
  • FILENAME_1: specifies the input file describing the descritization, following the same idea as for the ddm in the u-direction.
  • NAME_2: specifies the type of the finite elements space to allocate in the v-direction
  • FILENAME_2: specifies the input file describing the descritization, following the same idea as for the ddm in the v-direction.
  • NAME_3: specifies the type of the finite elements space to allocate in the w-direction
  • FILENAME_3: specifies the input file describing the descritization, following the same idea as for the ddm in the w-direction.

So, in this example, we are using a B-Spline discretization, with a tensor ddm.Next, let’s see what’s inside the two files:

> more input/parameters_ddm.nml
&P_DIMENSION
  DDM_DIM          = 3,
 /
&PARAMETERS
  PARTITIONER_NAME = 'TENSOR',
  N_PROCS          = 1, 1, 1,
 /

So if, you want to run the simulation in parallel, you’ll only need to change the values of N_PROCS to match the total number of processors you’re using when invoking MPI.

Now, let’s take a look to the other files:

> more input/parameters_bspline_u.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 16,
 DEGREES=          2,
 BC_MIN=           0,
 BC_MAX=           0,
 /

> more input/parameters_bspline_v.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 16,
 DEGREES=          2,
 BC_MIN=           0,
 BC_MAX=           0,
 /

> more input/parameters_bspline_w.nml
&P_DIMENSION
 P_DIM=          1,
 /
&PARAMETERS
 N_ELEMENTS= 16,
 DEGREES=          2,
 BC_MIN=           0,
 BC_MAX=           0,
 /

So as you can see, we are using a B-Spline finite elements discrezation on a grid of n_elements = 16 \times 16 \times 16 with a polynomial degrees (2,2,2).

Finally, the boundary conditions are specified by giving the code ID 0 which means that we are removing the interpolatory basis functions from the finite elements space.

Now let’s run our binary:

> bin/poisson_3d_ex01_01 input/parameters_sim_1d1d.nml input/solver_driver.nml
 >>>> context
 * assembly id          111
 * dimension            3
 * ddm name tensor
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 * discretization_2 bspline
 * discretization_3 bspline
 <<<<
 >>>> context
 * assembly id          111
 * dimension            3
 * ddm name tensor
 * ddm filename input/parameters_ddm.nml
 * discretization_1 bspline
 * discretization_2 bspline
 * discretization_3 bspline
 <<<<
 >>>>> parameters
 partitioner_id   :            1
 partitioner_name : TENSOR
 ddm_dim          :            3
 n_procs         :            1           1           1
 <<<<<
Assembly time =     4.412 seconds, on       4096 elements.
 * cg:  convegence after           5  iterations. Error    5.3970651295386319E-013
Assembly time =     0.464 seconds, on       4096 elements.
 >> norms :
   2.2625284449006593E-005   2.7788377734458527E-003

Visualization

We provide two ways for the visualization. The first one is based on the output object offred by FEMA which will be detailed in the next part. The result of the output export is a tecplot file that can be read by a visualization tool like VisIt or ParaView.

The second one is based on plot_field script available in DISCO. In the following example, we’ll directly export the field and the mapping into namelist files that we’ll process using Python. After executing our binary, we see the following new files:

phi.nml
mapping.nml

Now in order to visualize the solution stored in the file phi.nml, run:

> python $CLAPP_DIR/scripts/disco/plot_field.py phi.nml mapping.nml

In opposition to the 1d and 2d cases, you notice that the output is a vtk file that you need to open using your favorite tool.

The result is the following plot

alternate text

Numerical solution u_h

Dive into the code

It’s time to open the ex01_01.F90 file.

Initialize and Finalize subroutines

This part is exactly similar to the 1d case.

Declarations and polymorphic variables

This part is exactly similar to the 1d case.

Geometry and Mapping

In this example we are using a trilinear mapping between the default logical patch (a unit cube [0,1] \times [0,1] \times [0,1] ) and the cube \left( P_{ijk}, ~ i,j,k \in {1,2} \right)

! ... creates a trilinear mapping
P_111 = (/ 0.0, 0.0, 0.0 /)
P_211 = (/ 1.0, 0.0, 0.0 /)
P_121 = (/ 0.0, 1.0, 0.0 /)
P_221 = (/ 1.0, 1.0, 0.0 /)
P_112 = (/ 0.0, 0.0, 1.0 /)
P_212 = (/ 1.0, 0.0, 1.0 /)
P_122 = (/ 0.0, 1.0, 1.0 /)
P_222 = (/ 1.0, 1.0, 1.0 /)

call spl_mapping_trilinear( mapping, &
                          & P_111, P_121, P_211, P_221, &
                          & P_112, P_122, P_212, P_222)
call mapping % export('mapping.nml')
! ...

Maybe you noticed that we created our finite elements space without specifying the mapping. There are in fact three ways to use the mapping:

  1. create the space with the mapping.
  2. create the space then associate the mapping using the method set_mapping.
  3. create the space without a mapping, but call the assembler with the argument mapping.

In this example we choosed the 3 approach.

DDM

This part is exactly similar to the 1d case.

Linear Algebra

This part is exactly similar to the 1d case.

More about the assembler object

Now, we have our matrix and rhs object, but how to link them to the assembler object?

The assembler offers three methods to specify and associate objects like

  • matrices using the method set_matrix
  • rhs using the method set_rhs
  • fields using the method set_field

Please notice that these methods take an optional id that describes its position within the assembler data structure.

! ... associates the matrix, rhs and field to the assembler
call assembler % set_matrix(matrix, i_matrix=1)
call assembler % set_rhs(rhs, i_rhs=1)
call assembler % set_field(phi, i_field=1)
! ...

Now that we have associated all these objects, the only thing that remains is to specify the weak formulation and then assemble it.

Because we want to assemble two objects (matrix and rhs) we need to associate the assembler to two subroutines

  • one for the matrix and computes the integral between the trial and test basis functions v_j and v_i. This is done inside the subroutine build_matrix_stiffness.
  • one for the rhs and computes the integral between the function f and test basis function v_i. This is done inside the subroutine build_rhs.
! ... sets the weak formulation for both the matrix and rhs
assembler % ptr_matrix_contribution => build_matrix_stiffness
assembler % ptr_rhs_contribution    => build_rhs

call assembler % assemble(mapping=mapping)
! ...

Now let’s take a look to these subroutines. First, you’ll notice that build_matrix_stiffness is not defined in our example. In fact, FEMA provides some usefull weak formulations in the gallery package (mass, stiffness for example). You can find this subroutine in the file fema/fortran/src/gallery/gallery_3d.F90

! ..................................................
subroutine build_matrix_stiffness(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(jrk_rk) :: vi_0
  real(jrk_rk) :: vi_x
  real(jrk_rk) :: vi_y
  real(jrk_rk) :: vi_z
  real(jrk_rk) :: vj_0
  real(jrk_rk) :: vj_x
  real(jrk_rk) :: vj_y
  real(jrk_rk) :: vj_z
  real(jrk_rk) :: wvol

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % matrix_contribution = 0.0_jrk_rk
  ! ...

  ! ...
  do i_point = 1, n_points

    vi_0  = self % arr_vi(1,  1, i_point)
    vi_x  = self % arr_vi(1,  2, i_point)
    vi_y  = self % arr_vi(1,  3, i_point)
    vi_z  = self % arr_vi(1,  4, i_point)

    vj_0  = self % arr_vj(1,  1, i_point)
    vj_x  = self % arr_vj(1,  2, i_point)
    vj_y  = self % arr_vj(1,  3, i_point)
    vj_z  = self % arr_vj(1,  4, i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % matrix_contribution = self % matrix_contribution &
                             & + ( vi_x * vj_x + vi_y * vj_y + vi_z * vj_z ) * wvol
  end do
  ! ...

end subroutine build_matrix_stiffness
! ..................................................

You recover here our naming convention on the basis functions. You also notice how we compute the derivatives (vi_x for example). The aim of this subroutine (which is called for every couple (v_j, v_i) is to evaluate and cumulate the values of basis functions, their derivatives and some input data on the quadrature points, then update the attribut matrix_contribution (more details about this attribut will be given later).

Now, for the rhs, it is defined in our poisson file

! ..................................................
subroutine build_rhs(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: vi_0
  real(plf_rk) :: f_value
  real(plf_rk) :: wvol
  real(plf_rk) :: x
  real(plf_rk) :: y
  real(plf_rk) :: z
  real(kind=plf_rk), dimension(self % n_points) :: f_values

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  self % rhs_contribution = 0.0_plf_rk
  ! ...

  ! ...
  f_values = (k1**2 + k2**2 + k3**2) &
         & * sin(k1 * self % arr_x (1, :)) &
         & * sin(k2 * self % arr_x (2, :)) &
         & * sin(k3 * self % arr_x (3, :))
  ! ...

  ! ...
  do i_point = 1, n_points

    x     = self % arr_x (1, i_point)
    y     = self % arr_x (2, i_point)
    z     = self % arr_x (3, i_point)

    vi_0  = self % arr_vi(1,  1, i_point)

    f_value = f_values(i_point)

    wvol = self % weights(i_point) * self % jacobians(i_point)

    self % rhs_contribution = self % rhs_contribution + f_value * vi_0 * wvol
  end do
  ! ...

end subroutine build_rhs
! ..................................................

We recover the same accumulation principal that was introduced for the matrix. Here, the aim is to update the attribut rhs_contribution.

Linear Algebra Cont’

This part is exactly similar to the 1d case.

Computing the L^2 and H^1 norms

Following the same idea for the weak formulation, we will associate to the assembler object the following two subroutines

! ... associates fields and norms evaluation
assembler % ptr_fields_evaluation   => evaluate_fields
assembler % ptr_norms_evaluation    => evaluate_norms

call assembler % assemble(mapping=mapping)
! ...

The first subroutine evaluate_fields is used here to evaluate u_h (its first order derivative) on the quadrature points.

Note

This subroutine can be used if one deals with a nonlinear problem.

evaluate_fields is defined as follows

! ..................................................
subroutine evaluate_fields(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_trial_basis
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: vj_0
  real(plf_rk) :: vj_x
  real(plf_rk) :: vj_y
  real(plf_rk) :: vj_z

  ! ...
  do i_point = 1, self % n_points

    vj_0  = self % arr_vj(1, 1, i_point)
    vj_x  = self % arr_vj(1, 2, i_point)
    vj_y  = self % arr_vj(1, 3, i_point)
    vj_z  = self % arr_vj(1, 4, i_point)

    self % field_values(1, 1, :, i_point) = self % field_values(1, 1, :, i_point) &
                                  & + vj_0 * self % field_coeffs(:)

    self % field_values(1, 2, :, i_point) = self % field_values(1, 2, :, i_point) &
                                  & + vj_x * self % field_coeffs(:)

    self % field_values(1, 3, :, i_point) = self % field_values(1, 3, :, i_point) &
                                  & + vj_y * self % field_coeffs(:)

    self % field_values(1, 4, :, i_point) = self % field_values(1, 4, :, i_point) &
                                  & + vj_z * self % field_coeffs(:)
  end do
  ! ...

end subroutine evaluate_fields
! ..................................................

Finally the subroutine evaluate_norms is defined as

! ..................................................
subroutine evaluate_norms(self)
implicit none
  class(jrk_t_assembler_abstract), intent(inout) :: self
  ! local
  integer :: i_diagnostic
  integer :: i_field
  integer :: i_point
  integer :: n_points
  real(plf_rk) :: wvol
  real(plf_rk) :: norm_l2
  real(plf_rk) :: norm_h1
  real(plf_rk) :: f_0
  real(plf_rk) :: f_x
  real(plf_rk) :: f_y
  real(plf_rk) :: f_z
  real(plf_rk) :: field_0
  real(plf_rk) :: field_x
  real(plf_rk) :: field_y
  real(plf_rk) :: field_z
  real(kind=plf_rk), dimension(4, self % n_points) :: f_values

  ! ...
  n_points = self % n_points
  ! ...

  ! ...
  f_values(1, :) = sin(k1 * self % arr_x (1, :)) &
               & * sin(k2 * self % arr_x (2, :)) &
               & * sin(k3 * self % arr_x (3, :))
  f_values(2, :) = k1 * cos(k1 * self % arr_x (1, :)) &
                 &    * sin(k2 * self % arr_x (2, :)) &
                 &    * sin(k3 * self % arr_x (3, :))
  f_values(3, :) = k2 * sin(k1 * self % arr_x (1, :)) &
               &      * cos(k2 * self % arr_x (2, :)) &
               &      * sin(k3 * self % arr_x (3, :))
  f_values(4, :) = k3 * sin(k1 * self % arr_x (1, :)) &
               &      * sin(k2 * self % arr_x (2, :)) &
               &      * cos(k3 * self % arr_x (3, :))
  ! ...

  ! ...
  i_diagnostic = 0
  do i_field = 1, self % n_fields
    ! ...
    norm_l2 = 0.0_plf_rk
    norm_h1 = 0.0_plf_rk
    do i_point = 1, n_points

      f_0 = f_values(1, i_point)
      f_x = f_values(2, i_point)
      f_y = f_values(3, i_point)
      f_z = f_values(4, i_point)

      field_0 = self % field_values(1, 1, i_field, i_point)
      field_x = self % field_values(1, 2, i_field, i_point)
      field_y = self % field_values(1, 3, i_field, i_point)
      field_z = self % field_values(1, 4, i_field, i_point)

      wvol = self % weights(i_point) * self % jacobians(i_point)

      norm_l2 = norm_l2 +  (field_0 - f_0)**2 * wvol
      norm_h1 = norm_h1 + ((field_x - f_x)**2 + (field_y - f_y)**2 + (field_z - f_z)**2) * wvol
    end do
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_l2)
    ! ...

    ! ...
    i_diagnostic = i_diagnostic + 1
    call self % ptr_norms % add_value(i_diagnostic, norm_h1)
    ! ...
  end do
  ! ...

end subroutine evaluate_norms
! ..................................................

You notice that we are using the field_values array and then we update the pointer self % ptr_norms to the norms vector that we gave when creating the assembler object.

Output plot

The output object offred by FEMA allows to evaluate a field on the corresponding mesh, providing a given number of points per element.

! ..................................................
call output % create( export_name = "poisson_3d_ex01_01", &
                      & space     = trial_space , &
                      & n_fields  = 1           ,     &
                      & n_points_per_elements = 2 )


call output % set_field(phi, i_field=1)
call output % export()
! ..................................................

In this example, the export will generate the file poisson_3d_ex01_01_0.tp.

Memory deallocation

This part is exactly similar to the 1d case.