From 26d5cd3d2222859dbc77d50be7b30a653dab4111 Mon Sep 17 00:00:00 2001 From: arps Date: Mon, 14 Aug 2023 14:58:48 -0400 Subject: [PATCH 1/4] Initialized docker file --- Dockerfile | 31 +++++++++++++++++++++++++++++++ requirements.txt | 4 ++++ 2 files changed, 35 insertions(+) create mode 100644 Dockerfile create mode 100644 requirements.txt diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..c14f689 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,31 @@ +# Use the official Miniconda3 image as the base image +FROM continuumio/miniconda3 + +# Update package list and install necessary tools +RUN apt-get update && \ + apt-get install -y gfortran libblas-dev liblapack-dev + +# Set up an example Conda environment and install packages +RUN conda create -n myenv python=3.8 -y && \ + echo "source activate myenv" > ~/.bashrc + +# Set up the environment for running commands +SHELL ["/bin/bash", "-c", "source activate myenv"] + +# Install packages within the Conda environment +RUN conda install -c conda-forge numpy scipy matplotlib petsc4py -y + +# Set the environment variable +ENV EMCLAW_PATH /home/$USER/emclaw + +# Clone the repository and navigate to the riemann directory +RUN git clone https://github.com/MaxwellGEMS/emclaw.git $EMCLAW_PATH && \ + cd $EMCLAW_PATH/emclaw/riemann + +# Activate the Conda environment and build/install the package +RUN source activate myenv && \ + python setup.py build_ext -i && \ + conda develop -n myenv $EMCLAW_PATH + +# Set the default command to activate the Conda environment +CMD ["bash", "-i"] diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..2718b87 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +Flask==2.0.1 +requests>=2.25.1 +numpy==1.21.0 +pandas>=1.3.0 \ No newline at end of file From 7fd5e3fe8129901e3fb18ed26ca16721f062accc Mon Sep 17 00:00:00 2001 From: arps Date: Mon, 25 Sep 2023 15:26:00 -0400 Subject: [PATCH 2/4] Added comments to _convergence_sint --- examples/maxwell_vc_1d/_convergence_sint.py | 114 +++++++++++--------- 1 file changed, 62 insertions(+), 52 deletions(-) diff --git a/examples/maxwell_vc_1d/_convergence_sint.py b/examples/maxwell_vc_1d/_convergence_sint.py index 0b2fd56..8798270 100644 --- a/examples/maxwell_vc_1d/_convergence_sint.py +++ b/examples/maxwell_vc_1d/_convergence_sint.py @@ -1,14 +1,17 @@ +# Import necessary libraries import sys import os import numpy as np +# Adding paths to custom modules sys.path.append(os.path.realpath('../utils')) sys.path.append(os.path.realpath('../')) - +# Import custom modules from utils.materials import Material1D from utils.sources import Source1D +# Defining some constants and initialize material properties x_lower = 0. x_upper = 100. @@ -19,35 +22,40 @@ material.setup() material.delta_length = x_upper material.delta_corner = x_lower -material.delta_omega = 12.0*np.pi/200.0 +material.delta_omega = 12.0 * np.pi / 200.0 -source = Source1D(material,shape='off',wavelength=2.0) +# Create a source object +source = Source1D(material, shape='off', wavelength=2.0) source.setup() source.offset = 10.0 source.pulse_width = 2.0 -def grid_basic(x_lower,x_upper,mx,cfl): - dx = (x_upper-x_lower)/mx - dt = 0.9*cfl/(material.co*np.sqrt(1.0/(dx**2))) +# Define a function for grid calculations +def grid_basic(x_lower, x_upper, mx, cfl): + dx = (x_upper - x_lower) / mx + dt = 0.9 * cfl / (material.co * np.sqrt(1.0 / (dx ** 2))) tf = 100.0 - return dx,dt,tf + return dx, dt, tf -def dq_source(solver,state,dt): +# Define a function for source terms +def dq_source(solver, state, dt): """ Source terms for Maxwell's equations in one spatial dimension. Assume aux values are averages on the cell i """ - dq = -dt*state.q*state.aux[2:4,:]/state.aux[0:2,:] + dq = -dt * state.q * state.aux[2:4, :] / state.aux[0:2, :] return dq -def em1D(mx=1024,num_frames=5,cfl=1.0,outdir='./_output',before_step=True,debug=False,chi3=0.0,chi2=0.0,nl=False,psi=True,em=True,homogeneous=True): +# Define the main simulation function +def em1D(mx=1024, num_frames=5, cfl=1.0, outdir='./_output', before_step=True, debug=False, chi3=0.0, chi2=0.0, nl=False, psi=True, em=True, homogeneous=True): import clawpack.petclaw as pyclaw import petsc4py.PETSc as MPI + # Set material properties if not homogeneous if not homogeneous: if nl: material.chi3_e = chi3 @@ -56,46 +64,48 @@ def em1D(mx=1024,num_frames=5,cfl=1.0,outdir='./_output',before_step=True,debug= material.chi3_m = chi3 material.chi2_m = chi2 - if MPI.COMM_WORLD.rank==0: + if MPI.COMM_WORLD.rank == 0: material._outdir = outdir - source._outdir = outdir + source._outdir = outdir material._dump_to_latex() source._dump_to_latex() - num_eqn = 2 + num_eqn = 2 num_waves = 2 - num_aux = 4 + num_aux = 4 -# grid pre calculations and domain setup - dx,dt,tf = grid_basic(x_lower,x_upper,mx,cfl) - x = pyclaw.Dimension('x',x_lower,x_upper,mx) + # Grid pre-calculations and domain setup + dx, dt, tf = grid_basic(x_lower, x_upper, mx, cfl) + x = pyclaw.Dimension('x', x_lower, x_upper, mx) domain = pyclaw.Domain([x]) -# Solver settings - solver=pyclaw.SharpClawSolver1D() - solver.num_waves = num_waves - solver.num_eqn = num_eqn + # Solver settings + solver = pyclaw.SharpClawSolver1D() + solver.num_waves = num_waves + solver.num_eqn = num_eqn solver.weno_order = 5 solver.dt_variable = True - solver.dt_initial = dt/2.0 - solver.dt_max = dt - solver.max_steps = int(2*tf/dt) + solver.dt_initial = dt / 2.0 + solver.dt_max = dt + solver.max_steps = int(2 * tf / dt) -# Set the source + # Set the source if not psi: - print 'using dq_src' + # Python 2-style print statement + # print 'using dq_src' + print('using dq_src') solver.dq_src = dq_source -# Import Riemann and Tfluct solvers + # Import Riemann and Tfluct solvers if homogeneous: import maxwell_1d_rp else: import maxwell_1d_nl_rp as maxwell_1d_rp solver.tfluct_solver = False - solver.fwave = True - + solver.fwave = True + solver.rp = maxwell_1d_rp if solver.tfluct_solver: @@ -105,11 +115,11 @@ def em1D(mx=1024,num_frames=5,cfl=1.0,outdir='./_output',before_step=True,debug= import maxwell_1d_nl_tfluct as maxwell_1d_tfluct solver.tfluct = maxwell_1d_tfluct - - solver.cfl_max = cfl+0.05 + + solver.cfl_max = cfl + 0.05 solver.cfl_desired = cfl -# boundary conditions + # Boundary conditions solver.bc_lower[0] = pyclaw.BC.wall solver.bc_upper[0] = pyclaw.BC.wall @@ -118,47 +128,47 @@ def em1D(mx=1024,num_frames=5,cfl=1.0,outdir='./_output',before_step=True,debug= solver.reflect_index = [0] -# before step configure + # Before step configure if before_step: solver.call_before_step_each_stage = True solver.before_step = material.update_aux -# state setup - state = pyclaw.State(domain,num_eqn,num_aux) - + # State setup + state = pyclaw.State(domain, num_eqn, num_aux) + state.problem_data['chi2_e'] = material.chi2_e state.problem_data['chi3_e'] = material.chi3_e state.problem_data['chi2_m'] = material.chi2_m state.problem_data['chi3_m'] = material.chi3_m - state.problem_data['eo'] = material.eo - state.problem_data['mo'] = material.mo - state.problem_data['co'] = material.co - state.problem_data['zo'] = material.zo - state.problem_data['dx'] = state.grid.x.delta - state.problem_data['nl'] = nl - state.problem_data['psi'] = psi - - source._dx = state.grid.x.delta + state.problem_data['eo'] = material.eo + state.problem_data['mo'] = material.mo + state.problem_data['co'] = material.co + state.problem_data['zo'] = material.zo + state.problem_data['dx'] = state.grid.x.delta + state.problem_data['nl'] = nl + state.problem_data['psi'] = psi + + source._dx = state.grid.x.delta material._dx = state.grid.x.delta -# array initialization + # Array initialization source.init(state) material.init(state) - state.q = state.q*state.aux[0:2,:] - -# controller + state.q = state.q * state.aux[0:2, :] + + # Controller claw = pyclaw.Controller() claw.tfinal = tf claw.num_output_times = num_frames claw.solver = solver - claw.solution = pyclaw.Solution(state,domain) + claw.solution = pyclaw.Solution(state, domain) claw.outdir = outdir claw.write_aux_always = True return claw -if __name__=="__main__": +if __name__ == "__main__": import sys from clawpack.pyclaw.util import run_app_from_main - output = run_app_from_main(em1D) + output = run_app_from_main(em1D) \ No newline at end of file From df101768727eaae33cd1519066ac9b8cece5936d Mon Sep 17 00:00:00 2001 From: arps Date: Wed, 27 Sep 2023 14:15:54 -0400 Subject: [PATCH 3/4] Added comments to maxwell_1d.py --- examples/maxwell_vc_1d/maxwell_1d.py | 83 +++++++++++++++------------- 1 file changed, 45 insertions(+), 38 deletions(-) diff --git a/examples/maxwell_vc_1d/maxwell_1d.py b/examples/maxwell_vc_1d/maxwell_1d.py index 7dd8e16..cb6abbb 100755 --- a/examples/maxwell_vc_1d/maxwell_1d.py +++ b/examples/maxwell_vc_1d/maxwell_1d.py @@ -1,22 +1,27 @@ +# Import necessary modules import numpy as np from emclaw.utils.materials import Material1D from emclaw.utils.sources import Source1D from emclaw.utils import basics +# Create a 1D material with a homogeneous shape and set it up material = Material1D(shape='homogeneous') material.setup() -source = Source1D(material,shape='off',wavelength=2.0) +# Create a 1D source with shape 'off', wavelength 2.0, and set it up +source = Source1D(material, shape='off', wavelength=2.0) source.offset = 25.0 source.setup() x_lower = 0. x_upper = 100. -def em1D(mx = 1024, num_frames = 10, use_petsc = True, cfl = 1.0, conservative = True, - chi3 = 0.0, chi2 = 0.0, nl = False, psi = True, em = True, before_step = False, - debug = False, outdir = './_output', output_style = 1): +# Define a function em1D with various parameters +def em1D(mx=1024, num_frames=10, use_petsc=True, cfl=1.0, conservative=True, + chi3=0.0, chi2=0.0, nl=False, psi=True, em=True, before_step=False, + debug=False, outdir='./_output', output_style=1): + # Check if petsc4py is used, and import necessary modules if use_petsc: import clawpack.petclaw as pyclaw import petsc4py.PETSc as MPI @@ -30,39 +35,40 @@ def em1D(mx = 1024, num_frames = 10, use_petsc = True, cfl = 1.0, conservative = material.chi3_m = chi3 material.chi2_m = chi2 - if np.logical_and(use_petsc, MPI.COMM_WORLD.rank==0): - basics.set_outdirs(material, source, outdir = outdir, debug = debug) + # Set output directories based on MPI rank if using petsc4py + if np.logical_and(use_petsc, MPI.COMM_WORLD.rank == 0): + basics.set_outdirs(material, source, outdir=outdir, debug=debug) else: - basics.set_outdirs(material, source, outdir = outdir, debug = debug) + basics.set_outdirs(material, source, outdir=outdir, debug=debug) - num_eqn = 2 + num_eqn = 2 num_waves = 2 - num_aux = 4 + num_aux = 4 -# grid pre calculations and domain setup + # Grid precalculations and domain setup dx, dt, tf = basics.grid_basic([[x_lower, x_upper, mx]], cfl, material.co, source.v) - x = pyclaw.Dimension(x_lower,x_upper,mx,name='x') + x = pyclaw.Dimension(x_lower, x_upper, mx, name='x') domain = pyclaw.Domain([x]) -# Solver settings - solver=pyclaw.SharpClawSolver1D() - solver.num_waves = num_waves - solver.num_eqn = num_eqn + # Solver settings + solver = pyclaw.SharpClawSolver1D() + solver.num_waves = num_waves + solver.num_eqn = num_eqn solver.lim_type = 2 solver.dt_variable = True - solver.dt_initial = dt/2.0 - solver.dt_max = dt - solver.max_steps = int(2*tf/dt) + solver.dt_initial = dt / 2.0 + solver.dt_max = dt + solver.max_steps = int(2 * tf / dt) -# Import Riemann and Tfluct solvers + # Import Riemann and Tfluct solvers based on 'conservative' flag if conservative: from emclaw.riemann import maxwell_1d_rp else: from emclaw.riemann import maxwell_1d_nc_rp as maxwell_1d_rp solver.tfluct_solver = True - solver.fwave = True + solver.fwave = True solver.rp = maxwell_1d_rp @@ -74,10 +80,10 @@ def em1D(mx = 1024, num_frames = 10, use_petsc = True, cfl = 1.0, conservative = solver.tfluct = maxwell_1d_tfluct - solver.cfl_max = cfl+0.5 + solver.cfl_max = cfl + 0.5 solver.cfl_desired = cfl -# boundary conditions + # Boundary conditions solver.bc_lower[0] = pyclaw.BC.wall solver.bc_upper[0] = pyclaw.BC.wall @@ -86,43 +92,44 @@ def em1D(mx = 1024, num_frames = 10, use_petsc = True, cfl = 1.0, conservative = solver.reflect_index = [0] -# before step configure + # Before step configure if before_step: solver.call_before_step_each_stage = True solver.before_step = material.update_aux -# state setup - state = pyclaw.State(domain,num_eqn,num_aux) + # State setup + state = pyclaw.State(domain, num_eqn, num_aux) state.problem_data['chi2_e'] = material.chi2_e state.problem_data['chi3_e'] = material.chi3_e state.problem_data['chi2_m'] = material.chi2_m state.problem_data['chi3_m'] = material.chi3_m - state.problem_data['eo'] = material.eo - state.problem_data['mo'] = material.mo - state.problem_data['co'] = material.co - state.problem_data['zo'] = material.zo - state.problem_data['dx'] = state.grid.x.delta - state.problem_data['nl'] = nl - state.problem_data['psi'] = psi + state.problem_data['eo'] = material.eo + state.problem_data['mo'] = material.mo + state.problem_data['co'] = material.co + state.problem_data['zo'] = material.zo + state.problem_data['dx'] = state.grid.x.delta + state.problem_data['nl'] = nl + state.problem_data['psi'] = psi state.problem_data['conservative'] = conservative - source._dx = state.grid.x.delta + source._dx = state.grid.x.delta material._dx = state.grid.x.delta -# array initialization + # Array initialization source.init(state) material.init(state) + # Multiply state.q by state.aux[0:2,:] if 'conservative' is True if conservative: - state.q = state.q*state.aux[0:2,:] + state.q = state.q * state.aux[0:2, :] -# controller + # Controller claw = pyclaw.Controller() claw.tfinal = tf claw.num_output_times = num_frames claw.solver = solver - claw.solution = pyclaw.Solution(state,domain) + claw.solution = pyclaw.Solution(state, domain) claw.outdir = outdir claw.write_aux_always = True claw.output_style = output_style @@ -130,6 +137,6 @@ def em1D(mx = 1024, num_frames = 10, use_petsc = True, cfl = 1.0, conservative = return claw -if __name__=="__main__": +if __name__ == "__main__": from clawpack.pyclaw.util import run_app_from_main output = run_app_from_main(em1D) From f29674992c466edc0fd4a7525f8b133da1ee57c0 Mon Sep 17 00:00:00 2001 From: arps Date: Wed, 27 Sep 2023 14:19:20 -0400 Subject: [PATCH 4/4] Refactored code --- examples/maxwell_vc_1d/setplot.py | 49 +++++++++++++++---------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/examples/maxwell_vc_1d/setplot.py b/examples/maxwell_vc_1d/setplot.py index 38b8aa9..822d5a4 100755 --- a/examples/maxwell_vc_1d/setplot.py +++ b/examples/maxwell_vc_1d/setplot.py @@ -1,4 +1,3 @@ - """ Set up the plot figures, axes, and items to be done for each frame. @@ -9,7 +8,6 @@ import numpy as np - #-------------------------- def setplot(plotdata): #-------------------------- @@ -20,12 +18,12 @@ def setplot(plotdata): Output: a modified version of plotdata. """ - plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.clearfigures() # Clear any old figures, axes, items data # Figure for total field intensity and refractive index plotfigure = plotdata.new_plotfigure(name='I and n', figno=0) - # plot fields intensity + # Plot fields intensity plotaxes = plotfigure.new_plotaxes() plotaxes.axescmd = 'subplot(211)' plotaxes.xlimits = 'auto' @@ -36,10 +34,10 @@ def setplot(plotdata): plotitem.plot_var = ehfields plotitem.plotstyle = '-' plotitem.color = 'b' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':2,'markersize':1} - # plot refractive index + # Plot refractive index plotaxes = plotfigure.new_plotaxes() plotaxes.axescmd = 'subplot(212)' @@ -51,7 +49,7 @@ def setplot(plotdata): plotitem.plot_var = refind plotitem.plotstyle = '-' plotitem.color = 'g' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':2,'markersize':1} @@ -69,7 +67,7 @@ def setplot(plotdata): plotitem.plot_var = efield plotitem.plotstyle = '-' plotitem.color = 'b' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':1,'markersize':1} # Set up for axes in this figure: @@ -83,7 +81,7 @@ def setplot(plotdata): plotitem.plot_var = refind plotitem.plotstyle = '-' plotitem.color = 'g' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':1,'markersize':1} plotfigure = plotdata.new_plotfigure(name='H field', figno=2) @@ -97,7 +95,7 @@ def setplot(plotdata): plotitem.plot_var = hfield plotitem.plotstyle = '-' plotitem.color = 'r' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':1,'markersize':1} plotaxes = plotfigure.new_plotaxes() @@ -110,12 +108,12 @@ def setplot(plotdata): plotitem.plot_var = refind plotitem.plotstyle = '-' plotitem.color = 'g' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':1,'markersize':1} plotfigure = plotdata.new_plotfigure(name='u and n', figno=3) - # plot fields intensity + # Plot fields intensity plotaxes = plotfigure.new_plotaxes() plotaxes.axescmd = 'subplot(211)' plotaxes.xlimits = 'auto' @@ -126,10 +124,10 @@ def setplot(plotdata): plotitem.plot_var = energy plotitem.plotstyle = '-' plotitem.color = 'b' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':2,'markersize':1} - # plot refractive index + # Plot refractive index plotaxes = plotfigure.new_plotaxes() plotaxes.axescmd = 'subplot(212)' @@ -141,22 +139,22 @@ def setplot(plotdata): plotitem.plot_var = refind plotitem.plotstyle = '-' plotitem.color = 'g' - plotitem.show = True # show on plot? + plotitem.show = True # Show on plot? plotitem.kwargs = {'linewidth':2,'markersize':1} - # Parameters used only when creating html and/or latex hardcopy + # Parameters used only when creating HTML and/or LaTeX hardcopy # e.g., via visclaw.frametools.printframes: - plotdata.printfigs = True # print figures - plotdata.print_format = 'png' # file format - plotdata.print_framenos = 'all' # list of frames to print - plotdata.print_fignos = 'all' # list of figures to print - plotdata.html = True # create html files of plots? + plotdata.printfigs = True # Print figures + plotdata.print_format = 'png' # File format + plotdata.print_framenos = 'all' # List of frames to print + plotdata.print_fignos = 'all' # List of figures to print + plotdata.html = True # Create HTML files of plots? plotdata.html_homelink = '../README.html' - plotdata.latex = True # create latex file of plots? - plotdata.latex_figsperline = 3 # layout of plots - plotdata.latex_framesperline = 1 # layout of plots - plotdata.latex_makepdf = False # also run pdflatex? + plotdata.latex = True # Create LaTeX file of plots? + plotdata.latex_figsperline = 3 # Layout of plots + plotdata.latex_framesperline = 1 # Layout of plots + plotdata.latex_makepdf = False # Also run pdflatex? return plotdata @@ -180,4 +178,3 @@ def hfield(current_data): def ehfields(current_data): ehfi = np.sqrt((current_data.q[0,:]/current_data.aux[0,:])**2 + (current_data.q[1,:]/current_data.aux[1,:])**2) return ehfi -