In [7]:

    

%%file heateq.equelle

# Heat conduction with no boundary conditions or source terms.

# Physics that requires specification
k : Scalar = InputScalarWithDefault("k", 0.3) # Heat diffusion constant.

# @afr: time step strategy assumed given outside.
dt : Scalar = InputScalarWithDefault("dt", 0.5) # Time step length.

# u0 should be given initial values (in user input)
u0 : Collection Of Scalar On AllCells()
u0 = InputCollectionOfScalar("u0", AllCells())

# Trying to add support for Dirichlet boundaries
dirichlet_boundary : Collection Of Face Subset Of (BoundaryFaces())
dirichlet_boundary = InputDomainSubsetOf("dirichlet_boundary", BoundaryFaces())
dirichlet_val : Collection Of Scalar On dirichlet_boundary
dirichlet_val = InputCollectionOfScalar("dirichlet_val", dirichlet_boundary)

# Compute interior transmissibilities.
vol = |AllCells()|                                         # Deduced type:  Collection Of Scalar On AllCells()
interior_faces = InteriorFaces()                           # Deduced type:  Collection Of Face
first = FirstCell(interior_faces)                          # Deduced type:  Collection Of Cell On interior_faces
							   # Equivalent to: Collection Of Cell On InteriorFaces()
second = SecondCell(interior_faces)                        # Deduced type:  Same as for 'first'.
itrans : Collection Of Scalar On interior_faces = k * |interior_faces| / |Centroid(first) - Centroid(second)| 

# Compute boundary transmissibilities.
bf = BoundaryFaces()
bf_cells = IsEmpty(FirstCell(bf)) ? SecondCell(bf) : FirstCell(bf)
bf_sign = IsEmpty(FirstCell(bf)) ? (-1 Extend bf) : (1 Extend bf)
btrans = k * |bf| / |Centroid(bf) - Centroid(bf_cells)|

# Compute quantities needed for boundary conditions.
dir_sign = bf_sign On dirichlet_boundary

# Compute flux for interior faces.
computeInteriorFlux : Function(u : Collection Of Scalar On AllCells()) -> Collection Of Scalar On InteriorFaces()
computeInteriorFlux(u) = {
    -> -itrans * Gradient(u)
}

# Compute flux for boundary faces.
computeBoundaryFlux : Function(u : Collection Of Scalar On AllCells()) -> Collection Of Scalar On BoundaryFaces()
computeBoundaryFlux(u) = {
    # Compute flux at Dirichlet boundaries.
    u_dirbdycells = u On (bf_cells On dirichlet_boundary)
    dir_fluxes = (btrans On dirichlet_boundary) * dir_sign * (u_dirbdycells - dirichlet_val)
    # Extending with zero away from Dirichlet boundaries (i.e. assuming no-flow elsewhere).
    -> dir_fluxes Extend BoundaryFaces()
}

# Compute the residual for the heat equation.
computeResidual : Function(u : Collection Of Scalar On AllCells()) -> Collection Of Scalar On AllCells()
computeResidual(u) = {
    ifluxes = computeInteriorFlux(u)
    bfluxes = computeBoundaryFlux(u)
    # Extend both ifluxes and bfluxes to AllFaces() and add to get all fluxes.
    fluxes = (ifluxes Extend AllFaces()) + (bfluxes Extend AllFaces())
    # Deduced type: Collection Of Scalar On AllCells()
    residual = u - u0 + (dt / vol) * Divergence(fluxes)
    -> residual
}

# NewtonSolve takes a function (that should accept the primary variable as input) and the initial guess for the primary variable.
explicitu = u0 - computeResidual(u0)
u = NewtonSolve(computeResidual, u0)

Output("explicitu", explicitu)
Output("u", u)


    

    




Overwriting heateq.equelle

In [9]:

    

!/home/jse/projects/equelle/compiler/ec < heateq.equelle > heateq.cpp


    

In [10]:

    

!cat heateq.cpp


    

    




// This program was created by the Equelle compiler from SINTEF.

#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <cmath>
#include <array>

#include "EquelleRuntimeCPU.hpp"

void ensureRequirements(const EquelleRuntimeCPU& er);

int main(int argc, char** argv)
{
    // Get user parameters.
    Opm::parameter::ParameterGroup param(argc, argv, false);

    // Create the Equelle runtime.
    EquelleRuntimeCPU er(param);

    ensureRequirements(er);

    // ============= Generated code starts here ================

    const Scalar k = er.inputScalarWithDefault("k", double(0.3));
    const Scalar dt = er.inputScalarWithDefault("dt", double(0.5));
    const CollOfScalar u0 = er.inputCollectionOfScalar("u0", er.allCells());
    const CollOfFace dirichlet_boundary = er.inputDomainSubsetOf("dirichlet_boundary", er.boundaryFaces());
    const CollOfScalar dirichlet_val = er.inputCollectionOfScalar("dirichlet_val", dirichlet_boundary);
    const CollOfScalar vol = er.norm(er.allCells());
    const CollOfFace interior_faces = er.interiorFaces();
    const CollOfCell first = er.firstCell(interior_faces);
    const CollOfCell second = er.secondCell(interior_faces);
    const CollOfScalar itrans = (k * (er.norm(interior_faces) / er.norm((er.centroid(first) - er.centroid(second)))));
    const CollOfFace bf = er.boundaryFaces();
    const CollOfCell bf_cells = er.trinaryIf(er.isEmpty(er.firstCell(bf)), er.secondCell(bf), er.firstCell(bf));
    const CollOfScalar bf_sign = er.trinaryIf(er.isEmpty(er.firstCell(bf)), er.operatorExtend(-double(1), bf), er.operatorExtend(double(1), bf));
    const CollOfScalar btrans = (k * (er.norm(bf) / er.norm((er.centroid(bf) - er.centroid(bf_cells)))));
    const CollOfScalar dir_sign = er.operatorOn(bf_sign, er.boundaryFaces(), dirichlet_boundary);
    std::function<CollOfScalar(const CollOfScalar&)> computeInteriorFlux = [&](const CollOfScalar& u) -> CollOfScalar {
        return (-itrans * er.gradient(u));
    };
    std::function<CollOfScalar(const CollOfScalar&)> computeBoundaryFlux = [&](const CollOfScalar& u) -> CollOfScalar {
        const CollOfScalar u_dirbdycells = er.operatorOn(u, er.allCells(), er.operatorOn(bf_cells, er.boundaryFaces(), dirichlet_boundary));
        const CollOfScalar dir_fluxes = ((er.operatorOn(btrans, er.boundaryFaces(), dirichlet_boundary) * dir_sign) * (u_dirbdycells - dirichlet_val));
        return er.operatorExtend(dir_fluxes, dirichlet_boundary, er.boundaryFaces());
    };
    std::function<CollOfScalar(const CollOfScalar&)> computeResidual = [&](const CollOfScalar& u) -> CollOfScalar {
        const CollOfScalar ifluxes = computeInteriorFlux(u);
        const CollOfScalar bfluxes = computeBoundaryFlux(u);
        const CollOfScalar fluxes = (er.operatorExtend(ifluxes, er.interiorFaces(), er.allFaces()) + er.operatorExtend(bfluxes, er.boundaryFaces(), er.allFaces()));
        const CollOfScalar residual = ((u - u0) + ((dt / vol) * er.divergence(fluxes)));
        return residual;
    };
    const CollOfScalar explicitu = (u0 - computeResidual(u0));
    const CollOfScalar u = er.newtonSolve(computeResidual, u0);
    er.output("explicitu", explicitu);
    er.output("u", u);

    // ============= Generated code ends here ================

    return 0;
}

void ensureRequirements(const EquelleRuntimeCPU& er)
{
    (void)er;
}

In [19]:

    

!g++ -Wall -I/usr/include/eigen3 -I/home/jse/projects/equelle/examples/include  -std=c++11 heateq.cpp -o heateq -lopmautodiff -lopmcore -ldunecommon /home/jse/projects/equelle/examples/build/libequelle_rt.a


    

In [20]:

    

!./heateq


    

    




grid_dim not found. Using default value '2'.
ny not found. Using default value '1'.
dy not found. Using default value '1'.
nx not found. Using default value '6'.
dx not found. Using default value '1'.
linsolver not found. Using default value 'umfpack'.
output_to_file not found. Using default value 'false'.
verbose not found. Using default value '0'.
max_iter not found. Using default value '10'.
abs_res_tol not found. Using default value '1e-06'.
k not found. Using default value '0.3'.
dt not found. Using default value '0.5'.
u0_from_file not found. Using default value 'false'.
ERROR: The group '' does not contain an element named 'u0'.
terminate called after throwing an instance of 'Opm::parameter::ParameterGroup::NotFoundException'
  what():  std::exception
Aborted (core dumped)

In [ ]: