Solving Magnetic Field Using Floating Random Walk

[gt24]

Hi everyone,
I am working on the Fortran code to solve a magnetic field by using "Floating Random Walk" method.
I use a square domain for this case and show in the attachment. The external magnetic strength comes from left to right boundary.
Both top and bottom are assuming as insulation boundaries. I need some expert to help me to check if I the boundary conditions is a correct set up. Thanks!

       IF (X >= (10.0-ER)) THEN             ! BOUNDARY ON RIGHT
         SUM = SUM + 0.0
         GO TO 20
           ELSE
         END IF
        IF (Y >= (10.0-ER)) THEN             ! BOUNDARY ON TOP
         SUM = SUM + 0.0
         GO TO 20
         ELSE
        END IF

        IF ((X < ER) .AND. (Y < 10)) THEN    ! BOUNDARY ON LEFT
         SUM = SUM + 2.0                             ! {MAGNETIC STRENGTH}
        GO TO 20
         ELSE
        END IF

        IF ((X > ER) .AND. (X <=10.0-ER) .AND. (Y < ER)) THEN   ! BOUNDARY ON BOTTOM
         SUM = SUM + 0.0
        GO TO 20
         ELSE
        END IF
  20    CONTINUE

 

[BvU]

Hello and ,

(I am familiar with ftn, not with floating random walk)

Do we know what ER is ? Or, for that matter, where X and Y come from ?
(In short: what is this piece of code (program?, subroutine?) supposed to do and where does it fit in ?

Anyway, your code is hard to follow due to the go to statements -- perhaps using the else if form can make it more legible ?

 

[jbunniii]

In three out of the four cases, you're doing SUM = SUM + 0.0. What is the purpose of this?

 

[Baluncore]

Here is some Basic code for a floating random walk.

'=====================================================================
' Floating Random Walk; solves for potential at a 2D point within bounds
'=====================================================================
' Written in FreeBASIC; https://en.wikipedia.org/wiki/FreeBASIC
'
' Uses a Floating Random Walk to solve for potential in a 2D geometric
'   space, defined by combining circles of specified potential, with straight
'   lines that can have a linear potential gradient.
'
' Ref; Monte Carlo Methods for Electromagnetics - Matthew N.O. Sadiku - 2009 - CRC Press
'
' A random walk will not step across a closed boundary
' A wall opening wider than 2·Et, may allow a walk to escape to infinity
'
'=====================================================================
' Complex datatype = coordinate pair = point = 2D vector
Type Complex
    As Double x     ' real component
    As Double y     ' imaginary component
End Type

'---------------------------------------------------------------------
' overload some arithmetic operators for complex
Operator + ( Byval a As Complex, Byval b As Complex ) As Complex    ' addition
Return Type( a.x + b.x, a.y + b.y )
End Operator

Operator - ( Byval a As Complex, Byval b As Complex ) As Complex    ' subtraction
Return Type( a.x - b.x, a.y - b.y )
End Operator

Operator * ( Byval a As Complex, Byval b As Complex ) As Complex    ' vector rotation
Return Type( a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x )
End Operator

'---------------------------------------------------------------------
' distance between two points, the length of a line
Function pythagoras( Byref a As Complex, Byref b As Complex ) As Double
    Dim As Double dx = a.x - b.x, dy = a.y - b.y
    Return Sqr( dx*dx + dy*dy )
End Function

'=====================================================================
' define the objects used to build the geometry of potential bounds
'=====================================================================
Type a_circle   ' can be both an external and an internal boundary
    As Complex Ctr  ' centre
    As Double rad   ' radius
    As Double pot   ' potential
End Type

'---------------------------------------------------------------------
Type a_line ' with option for a potential gradient along the segment
    As Complex head     ' end 1
    As Complex tail     ' end 2
    As Complex cjuv     ' conjugate unit vector, to rotate points, precompute once
    As Double pot_head  ' potential at head end
    As Double pot_tail  ' potential at tail end
    As Double length    ' length of line, precompute once
End Type

'=====================================================================
' track the size of the arrays
Dim As Long m_circles = 100, m_lines = 100  ' maximum size of arrays
Dim As Long n_circles =   0, n_lines =   0  ' last valid entry in array

Dim As a_circle circles( 1 To m_circles )   ' index zero is invalid
Dim As a_line   lines( 1 To m_lines   )     '   use dynamic arrays ?

'=====================================================================
' Input some test data, load the arrays with geometry and potentials
'=====================================================================
' remember to check there is sufficient space in the array

' Matthew N. O. Sadiku - Elements of Electromagnetics-Oxford University Press (2018)
' page 244.  Chapter 6. Electrostatic Boundary-Value Problems

' example            exact
'   x        y         v
'  0.250    0.250     6.797
'  0.250    0.500    18.203
'  0.250    0.750    43.203
'
'  0.500    0.250     9.541
'  0.500    0.500    25.000
'  0.500    0.750    54.053
'
'              100V
'   y=1  ________________
'       |                |
'       |       3        |
'       |                |
'       |                |
'    0V | 4            2 | 0V
'       |                |
'       |                |
'       |       1        |
'       |________________|
'      o       0V       x=1

n_lines += 1    ' next index in line array
With lines( n_lines )
    .head = Type(  0,  0 )      ' 1, bottom
    .tail = Type(  1,  0 )
    .pot_head = 0
    .pot_tail = 0
End With

n_lines += 1
With lines( n_lines )
    .head = Type(  1,  0 )      ' 2, rhs
    .tail = Type(  1,  1 )
    .pot_head = 0
    .pot_tail = 0
End With

n_lines += 1
With lines( n_lines )
    .head = Type(  1, 1 )       ' 3, top
    .tail = Type(  0, 1 )
    .pot_head = 100
    .pot_tail = 100
End With

n_lines += 1
With lines( n_lines )
    .head = Type( 0, 1 )        ' 4, lhs
    .tail = Type( 0, 0 )
    .pot_head = 0
    .pot_tail = 0
End With

'---------------------------------------------------------------------

n_circles += 1      ' an outer bound circle, to detect an escape
With circles( n_circles )
    .Ctr = Type( 0, 0 )
    .rad = 5
    .pot = -1e20    ' make an escape very obvious
End With

'=====================================================================
' Precompute length and conjugated unit vector, those will be used
'   repeatedly to work out the relative position of step S to the line
Dim As Long i
For i = 1 To n_lines
    With lines( i )
        .length = pythagoras( .head, .tail )    ' save line length
        If .length < 1e-6 Then
            Print " Line##### has zero length. "; n_lines
            Sleep: Stop
        End If
        .cjuv = ( .tail - .head )   ' translate line to place head at origin
        .cjuv.x /= + .length        ' divide by length scales to a unit vector
        .cjuv.y /= - .length        ' negate imaginary to conjugate
    End With
Next i

'---------------------------------------------------------------------
' avoid an empty data set
If ( n_circles + n_lines ) = 0 Then
    Print " Geometry data is missing. "
    Sleep : Stop
End If

'=====================================================================
' Geometry is specified, do the random walks
'=====================================================================
Const As Double BIG = 1e20      ' bigger than real data
Const As Double ET = 0.001      ' exit tolerance
Const As Long n_walks = 2000    ' how many walks to take

Const As Double Two_Pi = 8 * Atn( 1 )
Randomize( Timer * 1e6 )    ' seed the generator

Dim As Long walk, steps ' count the walks and steps
Dim As Long   closest_circle, closest_line  ' remember for walk's end
Dim As Double close_circle_r, close_line_r
Dim As Double sum, range, theta, r  ' sum of potentials, range to object
Dim As Complex So, S, T ' the start point, step, and the transformed position of S

So = Type( 0.5, 0.5 )   ' the point to analyse, origin of the random walks

'---------------------------------------------------------------------
steps = 0
sum = 0
For walk = 1 To n_walks
   
    S = So  ' start at the same point for each successive walk
    Do  ' for every step
        r = BIG     ' r will reduce to closest boundary
       
        '-------------------------------------------------------------
        ' find closest circle to S
        closest_circle = 0
        close_circle_r = BIG
        For i = 1 To n_circles
            With circles( i )
                range = pythagoras( S, .Ctr )   ' distance from centre
               
                If range > .rad Then
                    range = range - .rad    ' S is outside the circle
                Else
                    range = .rad - range    ' S is inside the circle
                End If
               
                If r > range Then r = range ' the minimum
               
                If close_circle_r > range Then
                    close_circle_r = range  ' remember range and
                    closest_circle = i      '   index of closest circle
                End If
               
            End With
        Next i
       
        '-------------------------------------------------------------
        ' find closest line to S
        closest_line = 0
        close_line_r = BIG
        For i = 1 To n_lines
            With lines( i )
               
                ' consider translating .head to origin, rotate to put .tail on +x axis
                T = ( S - .head ) * .cjuv   ' then S would be mapped to T
               
                If T.x < 0 Then             ' S is before head
                    range = pythagoras( S, .head )
                Else
                    If T.x > .length Then   ' S is beyond tail
                        range = pythagoras( S, .tail )
                    Else                    ' S is to side of line
                        range = Abs( T.y )
                    End If
                End If
               
                If r > range Then r = range ' the minimum
               
                If close_line_r > range Then
                    close_line_r = range    ' remember range and
                    closest_line = i        '   index of closest line
                End If
               
            End With
        Next i
       
        '-------------------------------------------------------------
        ' r is now the range S from the closest bound
        If r < ET Then Exit Do  ' exit terminal condition, so this walk ends
       
        '-------------------------------------------------------------
        ' not yet terminal, so take a step of length r in a random direction
        theta = Rnd * TWO_PI    ' random, 0 to 360 deg
        S.x += r * Cos( theta ) ' polar to rectangular
        S.y += r * Sin( theta ) '   advance to next step
       
        steps += 1  ' count steps over all walks to average
    Loop    ' go check this new step position
   
    '-----------------------------------------------------------------
    ' walk terminating, so identify closest boundary and evaluate potential
    If close_circle_r < close_line_r Then
        ' it was a circle boundary
        sum += circles( closest_circle ).pot
    Else
        ' it was a line boundary
        With lines( closest_line )
            T = ( S - .head ) * .cjuv   ' S will be mapped to T
            If T.x < 0 Then             ' S is before end 1
                sum += .pot_head
            Else
                If T.x > .length Then   ' S is beyond end 2
                    sum += .pot_tail
                Else    ' S is broadside to the line, interpolate potential
                    sum += ( .pot_tail * T.x + .pot_head * ( .length - T.x ) ) / .length
                End If
            End If
        End With
    End If
   
    '-----------------------------------------------------------------
Next walk   ' take another walk

'---------------------------------------------------------------------
Print
Print Using " x =##.##  x =##.##"; So.x; So.y
Print
Print Using " Potential =###.###"; sum / n_walks
Print
Print Using " Average steps per walk =####.###"; steps / n_walks
Print
Print " Done."

Sleep

'=====================================================================

 

[Baluncore]

In three out of the four cases, you're doing SUM = SUM + 0.0. What is the purpose of this?

When a walk ends on a boundary, it accumulates the potential of that boundary cell and uses it in the average for the walks from that start cell. If the top and bottom bounds are insulators, and the right is zero, then the left hand side only will have a non-zero potential. The accumulation of a zero, signals a zero potential edge, or an insulator. In this case, that occurs in three out of the four possible terminal conditions for a random walk originating from the point being solved.

The boundary cells of the rectangular potential array will contain the initial boundary conditions. There must be some way to flag an insulator that will have no influence on potential. The body of the array will then be filled with statistically computed potential values.

Any walk that ends at an insulated boundary cell must be ignored completely. Do not accumulate any potential, do not count it as a valid walk. It never happened.

 

[jbunniii]

When a walk ends on a boundary, it accumulates the potential of that boundary cell and uses it in the average for the walks from that start cell. If the top and bottom bounds are insulators, and the right is zero, then the left hand side only will have a non-zero potential. The accumulation of a zero, signals a zero potential edge, or an insulator. In this case, that occurs in three out of the four possible terminal conditions for a random walk originating from the point being solved.

The boundary cells of the rectangular potential array will contain the initial boundary conditions. There must be some way to flag an insulator that will have no influence on potential. The body of the array will then be filled with statistically computed potential values.

Any walk that ends at an insulated boundary cell must be ignored completely. Do not accumulate any potential, do not count it as a valid walk. It never happened.

I'm familiar with random walks, including random walks with boundaries.

But why is it necessarily to explicitly include the lines SUM = SUM + 0.0? Unless I'm missing some floating point subtlety (necessarily non-IEEE compliant!), these lines have no effect and the function would do exactly the same thing if they were omitted. And, especially because there are no comments explaining these lines, anyone reviewing or maintaining the code is likely to think that they are mistakes.

 

[Baluncore]

But why is it necessarily to explicitly include the lines SUM = SUM + 0.0?

It is quite clear to me what the posted code fragment is supposed to be doing. The presence of those lines as position holders, makes the structure obvious. Removal of the lines would have no effect.

There are several problems with the code structure. While not an experienced programmer, gt24 was willing to have a go. gt24 must have known those lines would have no effect during development. I expect gt24 intended to replace or modify them later. In hindsight it would have been better to insert a comment, rather than to use a dummy line that would draw unproductive criticism. Neophyte programmers need assistance, not criticism.