Hey,
i have solved it in my ReLU NN in the following way:
class ReLU_NN(tr.nn.Module):
    '''
    Class for a ReLU-NN with variable size.
    Input:
        -nn_list = [input dim., first hidden layer size,...,
                    last hidden layer size, output dim.]

    '''
    def __init__(self, nn_list):
        super(ReLU_NN, self).__init__()
        self.nn_list = nn_list
        self.hidden  = tr.nn.ModuleList()
        for i in range(len(nn_list)-1):
            self.hidden.append(tr.nn.Linear(nn_list[i], nn_list[i+1]).double())

    def forward(self, x):
        # forward pass through the network
        for layer in self.hidden[:-1]:
            x = tr.nn.functional.relu(layer(x))

        # the last layer is a linear layer without activation function
        output = self.hidden[-1](x)

        return output

    def add_layers(self, nn_add_list):
        '''
        Adds some layer between the last hidden layer and the output layer
        Input:
            -nn_add_list: a list which length defines the number of layer to be
                          added and each entry defines the number of neurons
                          e.g. [num1,...,numL]
                          Then self.hidden and self.nn_list will be adjusted
                          in the following way:
                          self.hidden = [old layer,...,old layer,new layer 1,...,
                                         new layer L, new output layer]
                          self.nn_list = [input-dim, old-layer-size,...,
                                          old-layer-size, num1,...,numL, output-
                                          dim]
        '''
        # adjust nn_list
        length = len(nn_add_list)
        temp   = self.nn_list[-1]
        self.nn_list.extend(nn_add_list)
        for i in range(length, 0, -1):
            self.nn_list[-i-1] = self.nn_list[-i]
        self.nn_list[-1] = temp

        # adjust hidden
        for i in range(length+2, 2, -1):
            self.hidden.insert(len(self.hidden)-1,\
                    tr.nn.Linear(self.nn_list[-i], self.nn_list[-i+1]).double())
        
        # adjust last hidden layer
        self.hidden[-1] = tr.nn.Linear(self.nn_list[-2], self.nn_list[-1]).double()

If have also implemented a method for adding neurons to the layers, but this is a little bit trickier. If someone wants it i can show the code snipped here. But the real problem is, that you have to change dynamically the parameter_groups within the optimizer and this is, as far as i know not possible. All you can do is to redefine the optimizer with the new model.parameters() (of course with the preoptimized weights and biases) and restart the optimization prozess.