Compile FORTRAN to WebAssembly and Solve Electromagnetic Fields in Web Browsers

Electric Field of a Microstrip Transmission Line

Pictured: Electric field lines of a microstrip transmission line over a ground plane. Image from Grounding and Shielding: Circuits and Interference by Ralph Morrison.

Introduction

In this article, a method of cross-compiling FORTRAN code to WebAssembly with LLVM and flang is described.

Since WebAssebly uses 32-bit pointers, but 32-bit platforms are not supported by flang, currently, only numerical computation is supported, it’s not possible to use the standard library in FORTRAN code, such as I/O (see Appendix 0: See Also for an alternative solution). However, in spite of this limitation, it already enabled a non-trivial practical application: cross-compiling the TNT-MMTL 2D electromagnetic field solver to WebAssembly, allowing electronics designers to solve microwave transmission lines using the Boundary Element Method (BEM) without leaving the Web browser.

It’s worth noting that most FORTRAN 77 subroutines used in TNT-MMTL came directly from the original LINPACK from 1978, written by Jack Dongarra and Cleve Moler - they have been successfully complied to WebAssembly to solve systems of equations in Web browsers without any modification. I’d call it a technical marvell.

Part 1: Toolchain

WebAssembly

WebAssembly is a technology that allows the execution portable binary code in a Web browser within a virtual machine sandbox via a hypothetical CPU and instruction set. For old-timers, it resembles a Pascal p-code machine. Recent versions of the LLVM compiler introduced a WebAssembly backend, allowing one to cross-complie code directly to any programming language front-end supported by LLVM, as if it’s just another CPU.

flang

There exists several different out-of-tree projects and forks to implement a FORTRAN front-end on top of LLVM, confusingly, all of them are named flang. Since LLVM 11 (2020), a flavor of flang has been merged into the upstream to become the official implementation. In this article, flang refers to the upstreamed version, it’s currently known as flang-new in the command-line.

Early flang-new versions were only capable of syntax checking, it was unable to generate any object code. Thus, the latest version of LLVM and flang must be used. At the time of writing, it’s LLVM 15 (2023). If your operating system doesn’t provide the latest version via the package manager, it’s necessary to build LLVM from source.

On Debian Sid, to install LLVM 15 with flang front-end and the LLVM linker lld, enter:

sudo apt install llvm-15
sudo apt install flang-new-15
sudo apt install lld-15

Note that for anything but the most trivial program, it’s also necessary to install clang runtime and C runtime libraries, as flang generates code that calls C, such as memset(). On Debian, these packages are installed by default as dependencies, but if a self-built LLVM is not configured correctly, there can be mysterious undefined references if these libraries are missing. See Appendix 1: libclang_rt and C runtime for a solution.

By default, these Debian binaries are installed with a version suffix to avoid clashing with other versions. To avoid these superfluous suffixes, add the underlying /usr/lib/llvm-15/bin to your search path.

export PATH="/usr/lib/llvm-15/bin:${PATH}"

First FORTRAN program - Return an Integer

To ensure that the toolchain is functioning correctly, we can now try cross-compiling our first FORTRAN program, test.f90.

function test() result(x) bind(c, name='test')
    integer :: x
    x = 42
end function test

This program defines a function test() that returns an integer 42. It also exports this function as a visible C function, for our purpose, it makes the function callable from JavaScript.

Build

Since the flang compiler currently doesn’t have any supporting code for the WebAssembly target, the cross-compiling must be performed manually in 3 steps:

  1. Use flang-new to translate FORTRAN code .F or .f90 to LLVM IR code .ll.

    The command-line switches -emit-llvm -S ask LLVM to generate assembly (-S) in LLVM IR code. Argument --target=i386-unknown-linux ensures that flang is targetting 32-bit platforms - though it currently doesn’t really make a difference, as support for 32-bit platform is incomplete.

    flang-new-15 --target=i386-unknown-linux -emit-llvm -S test.f90 -o test.ll

  2. Use llc to translate LLVM IR code .ll to 32-bit WebAssembly object code .o.

    llc --march=wasm32 -filetype=obj test.ll -o test.o

  3. Use wasm-ld to link 32-bit WebAssembly object code .o to .wasm executable.

    The command-line switch --no-entry allows linking to an executable without any entry points, like main(), and --export-all forces the linker to export all functions as external symbols.

    wasm-ld --no-entry --export-all test.o -o test.wasm

Run

To run the WebAssembly code in a Web browser, we create a test webpage test.html with the following content.

<!DOCTYPE html>
<html>
<body>
<script type="module">
    async function init() {
        const { instance } = await WebAssembly.instantiateStreaming(
            fetch("./test.wasm")
        );

        let result = instance.exports.test();
        document.write(result)
    }

    init();
</script>

</body>
</html>

As seen, it contains a few lines of JavaScript code that imports the WebAssembly executable ./test.wasm and calls the exported function test() from JavaScript, and writes the return value to HTML DOM.

Due to security reasons, web pages are not allowed to access arbitrary local files. This web page must be provided via a web server. Python has a built-in Web server which is often handy for this purpose.

python -m http.server

If the toolchain is functioning correctly, you should be able to see 42 after pointing a modern Web browser to the page. Here’s an example with Firefox.

Firefox, showing a webpage that says &ldquo;42&rdquo;

Second FORTRAN program - Leibniz’s formula for π

Next, let’s try a numerical computation problem - calculating a rough approximation of π via Leibniz’s formula by evaluating the following simple infinite series:

1−13+15−17+19−…=π4

This can be implemented by the following FORTRAN code, pi.f90.

module m_pi
    implicit none
    integer, parameter :: dp = selected_real_kind(15, 307)

contains

    ! calculate pi using Gregory-Leibniz's formula
    ! https://craftofcoding.wordpress.com/2020/04/09/calculating-π-with-gregory-leibniz-ii-fortran/
    function pi_leibniz(terms) result(x) bind(c, name='pi_leibniz')
       integer, value :: terms
       integer :: i, sign
       real (kind=dp) :: x
    
       x = 1.0_dp
    
       sign = -1
       do i = 1, terms
           x = x + (1.0_dp / ((i * 2.0_dp) + 1.0_dp)) * sign
           sign = sign * -1
       end do
    
       x = x * 4.0_dp
    end function pi_leibniz

end module m_pi

Build

flang-new -O2 --target=i386-unknown-linux -emit-llvm -S pi.f90
llc --march=wasm32 -filetype=obj pi.ll
wasm-ld  --no-entry --export-all -o pi.wasm pi.o

Run

Save the following code as pi.html and host it via a web server.

<!DOCTYPE html>
<html>
<body>
<script type="module">
    async function init() {
        const { instance } = await WebAssembly.instantiateStreaming(
            fetch("./pi.wasm")
        );

        let result = instance.exports.pi_leibniz(100000000);
        document.write(result)
    }

    init();
</script>

</body>
</html>

After opening this page, one should see the following result. Not all digits are correct, since Leibniz’s formula converges extremely slowly. This formula also has the unusual mathematical property that when this series is truncated, not all digits are wrong beyond a point. Instead, correct digits are interleaved with incorrect digits. This is to be expected and it’s not a bug.

Firefox, showing a webpage that shows the result of π

Part 2: Porting Field Solver TNT-MMTL to the Web

Motivation

TNT-MMTL is a 2D electromagnetic field solver that solves transmission lines with an arbitrary number of dielectric layers using the Boundary Element Method (a.k.a. Method of Moments) according to the laws of electrostatics.

It was developed by Mayo Clinic’s Special Purpose Processor Development Group (SPPDG) from the 1980s to the 1990s. In the early 2000s, it was released as free software under the GPLv2+ license. Development was then discontinued, and the project fell into obscurity and is mostly forgotten.

A field solver like TNT-MMTL can be used to find the characteristic impedance of a transmission line on a printed circuit board, which is of crucial importance in radio and high-speed digital electronics design. Considering that the industry standard solver, Si8000/Si9000 by Polar Instruments, is not just proprietary but also with a licensing fee of many thousand dollars, bringing TNT-MMTL to the Web would be an interesting idea. Although it doesn’t have as many features (or even generate the correct solution for some structures), still, for simple cases, the results are comparable within a percent or two.

Architecture

TNT-MMTL consists of two parts, the graphical user interface TNT, written in Tcl/Tk, and the Boundary Element Method field solver engine - BEM-MMTL, or just BEM. The BEM engine is written in C++, but it uses a FORTRAN library, Naval Surface Warfare Center Mathematical Library (NSWC) as its linear algebra kernel, it’s an old version written in FORTRAN 77. The code involved is truly ancient, most used subroutines came directly from the original LINPACK release from 1978, written by Jack Dongarra and Cleve Moler.

While it shouldn’t be too difficult to replace NSWC/LINPACK with another linear algebra library, or using f2c to translate them to C code, a faithful cross-complie is much preferred.

This can be achieved in three steps:

  1. Cross-compile all FORTRAN code to WebAssembly object code via flang, and generate a nswc.a library archive.
  2. Cross-compile all C++ code to WebAssembly object code via emscripten.
  3. Link both C++ and FORTRAN via emscripten.

Bypassing autotools

BEM uses GNU autotools as its build system, it was a reasonable choice - the project was developed in the 1980s, you can find mentions about compatibility problems on DEC Ultrix in its source code. But this creates several problems today for our hack.

First, the autotools files have not been updated since the early 2000s, and it’s seriously broken today (it was already broken at the time of release). It needs a big overhaul before it can be restored to work reliably. Next, autotools would still be incompatible with our highly non-standard build procedures anyway.

As the size of the codebase is small, it’s practical to bypass autotools and building all files from tree by invoking the compiler directly.

Cross-compiling FORTRAN code

Outside the source code tree ./src, I created a working directory called ./build-fortran, and wrote the following Makefile.

FLANG=flang-new
FFLAGS=-O2
LLC=llc
LD=wasm-ld
AR=llvm-ar

SRC_DIR ?= ../src

# List of FORTRAN files.
.ftrn_src_noprefix = \
        fft.F corth.F ceigv.F cmtms.F mtms.F cmslv1.F mslv.F dmslv.F \
        dcmslv.F sgefa.F sgeco.F dgefa.F sgesl.F dgesl.F fmin.F daxpy.F \
        ddot.F saxpy.F sdot.F idamax.F dscal.F sasum.F sscal.F isamax.F \
        dcfact.F dcsol.F dcminv.F dgeco.F dgedi.F cbal.F sgedi.F cgefa.F \
        cgesl.F comqr2.F cbabk2.F sfft.F spmpar.F cdivid.F dswap.F dasum.F \
        cgedi.F sswap.F caxpy.F cdotc.F icamax.F cscal.F ipmpar.F cswap.F \
        dpmpar.F

# Add $SRC_DIR prefix to all files to allow out-of-tree build
.ftrn_src = $(addprefix $(SRC_DIR), $(.ftrn_src_noprefix))

# Replace .F to .o to get output object code file
.ftrn_obj = $(.ftrn_src_noprefix:.F=.o)

# All object code files are archived as nswc.a library
nswc.a: $(.ftrn_obj)
	$(AR) rcs nswc.a $(.ftrn_obj)

# Invoke flang to build FORTRAN to WebAssembly
%.o: $(SRC_DIR)/%.F
	$(FLANG) $(FFLAGS) --target=i386-unknown-linux -emit-llvm -S $(SRC_DIR)/$*.F
	$(LLC) --march=wasm32 -filetype=obj $*.ll

clean:
	rm -rf *.ll *.o *.wasm *.a

After typing make, all FORTRAN files are cross-compiled to WebAssembly.

Cross-compiling C++ code

C++ code is cross-compiled in another directory build-cpp with another similar Makefile, and using emscripten (em++) as the compiler.

CXX=em++
# CXXFLAGS=
# omitted, list of CXXFLAGS extracted from autotools output

SRC_DIR ?= ../src

# .cpp_src_noprefix =
# omitted, list of C++ files

.cpp_src = $(addprefix $(SRC_DIR), $(.cpp_src_noprefix))

.cpp_obj = $(.cpp_src_noprefix:.cpp=.o)

bem.js: $(.cpp_obj)
	$(CXX) $(.cpp_obj) nswc.a -o bem.js --preload-file coplanar.xsctn -sEXPORTED_RUNTIME_METHODS=callMain,FS

%.o: $(SRC_DIR)/%.cpp
	$(CXX) -c $(CXXFLAGS) $(SRC_DIR)/$*.cpp -o $*.o

clean:
	rm -rf *.ll *.o *.wasm *.a

When bem.js is built, the FORTRAN library nswc.a is linked against it along with all C++ object files. Furthermore, emscripten functions are private by default, thus two functions callMain and FS are exported manually, allowing to call the main() function or manipulating the filesystem in JavaScript.

Finally, a test file coplanar.xsctn, the simulation file of a coplanar waveguide, is “preloaded” into the virtual filesystem for testing. I also slightly modified the BEM source code to print the simulation result to stdout.

Running in Web browser

Saving the following file to bem.html and start a Web server.

<script type="text/javascript">
function print(stdout)
{
  const escapeTrick = document.createElement('textarea');
  escapeTrick.textContent = stdout;
  const escaped = escapeTrick.innerHTML;
  
  document.write(escaped);
  document.write("<br>");
}

var Module = { 
    // Don't run main() on page load
    noInitialRun: true,    
    onRuntimeInitialized: () => {
        Module.callMain(["coplanar", "50", "50"])
    },  
    print: print,
    printErr: print,
};
</script>

<script src="bem.js"></script>

When opening bem.html from a Web browser, the following result can be seen:

Firefox, showing a webpage that shows the run log of BEM:filename: coplanarcntr_seg: 50 pln_seg: 50 coupling: 0.0254risetime: 2.5e-11 conductivity: 5.8e+07 frequency: 1e+08half_min_dim: 2.54e-06 grnd_planes: 1 top_grnd_thck: 0bot_grnd_thck: 1 num_sig: 1 num_grounds: 2units: 032719 MHz is the minimum frequency for the surface current assumptions for this cross section194 elements and 392 nodes were generatedlargest matrix to be inverted is 392 X 392Calculate LHS (assemble) matrix in free spacecalculate RHS (load) matrix for conductor 1Solve system of equationsIntegrate charge densityCalculate LHS (assemble) matrix in dielectriccalculate RHS (load) matrix for conductor 1Solve system of equationsIntegrate charge density Firefox, showing a webpage that shows the result of BEM:Mutual and Self Inductance:L(Active Signal , Passive Signal) Henrys/MeterL( ::condR0 , ::condR0 )= 3.2094513e-07Characteristic Impedance (Ohms):For Signal Line ::condR0= 30.7689Effective Dielectric Constant:For Signal Line ::condR0= 9.79222Propagation Velocity (meters/second):For Signal Line ::condR0= 9.5803224e+07Propagation Delay (seconds/meter):For Signal Line ::condR0= 1.0438062e-08Rdc:Rdc(Active Signal , Passive Signal) Ohms/MeterRdc( ::condR0 , ::condR0 )= 5.7407524e+01

The electromagnetic field solver MMTL has successfully obtained the characteristic impedance of the coplanar waveguide, the FORTRAN linear algebra code also ran flawlessly without any modification.

The full source code, including WebAssembly build instructions, can be found here:

Making it Practical

Based on this result, I’m now working on a practical version of this application. The Tcl/Tk GUI “TNT” will be reimplemented in JavaScript, as “WebTNT”. The first prototype is nearly-complete and it’s expected to be released within a week.

Here’s a quick demo (the GIF animation may take a while to load).

Firefox, showing the use of the graphical user interface to solve a microstrip and a coplanar waveguide transmission line interactively.

Furthermore, TNT-MMTL requires users to model each layer of the PCB manually. A further step forward is writing a GUI that works at a level higher than TNT, instead of modeling each layer of the board, it would generate the stackup automatically based on the type of transmission lines and its parameters, in Si8000/Si9000 style. My hope is to create the ultimate PCB impedance calculator on the Web.

Appendix 0: See Also

This article is a survey of all possible pathways from FORTRAN to WebAssembly as of 2020. The author developed a creative solution of compiling FORTRAN to GCC GIMPLE code and translating GIMPLE to LLVM IR. This entirely bypasses the 32-bit support problem of flang, allowing one to run full FORTRAN programs in WebAssembly, including standard libraries.

However, this method is not used as the incomplete flang solution already achieved my goals.

In Japanese. As far as I know, this was the only complete tutorial of FORTRAN to WebAssembly via flang, which was also used as the basis of my article.

Appendix 1: libclang_rt and C runtime

The front-end flang reuses clang code for some functionalities, this includes library libclang_rt (also known as compiler-rt), the built-in runtime library of clang. The C standard library is also needed for anything but the most trivial code (for example, memset() calls are often generated).

If you get mysterious undefined references during linking, it means LLVM was not built correctly. As a workaround, pre-built libraries can be downloaded from wasi-sdk project’s releases.

To install libclang_rt, download libclang_rt.builtins-wasm32-wasi-19.0.tar.gz, extract it to obtain lib/wasi/libclang_rt.builtins-wasm32.a.

Then pass

-Llib/wasi/libclang_rt.builtins-wasm32.a \
-l:libclang_rt.builtins-wasm32.a

when wasm32-ld is invoked. Similarly, missing libc libraries can be obtain by downloading wasi-sysroot-19.0.tar.gz and linking against .a files manually.